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SECTION  1 
INTRODUCTION 


The  WEDCOM  code  (Weapon  Effects  on  D-region  Communications)  is  a  digital 
computer  program  that  provides  calculations  of  ELF,  VLF,  and  LF  electric  and  magnetic 
field  strengths  in  ambient  and  nuclear  disturbed  environments  (Reference  1).  The 
code  is  intended  for  use  when  a  relatively  detailed  analysis  of  the  propagation  be¬ 
tween  two  terminals  is  required  in  nuclear  disturbed  environments.  The  computa¬ 
tional  models  (particularly  the  nuclear  environment  models)  are  intermediate  between 
first  principle  and  simplified  simulation  models.  Most  physical  quantities  are  ex¬ 
plicitly  modeled.  This  report  describes  work  on  code  structure  and  propagation  models 
in  order  to  reduce  link  evaluation  times. 

1.1  BACKGROUND  INFORMATION  AND  HISTORY 

OF  THE  WEDCOM  CODE 

Propagation  at  frequencies  in  the  VLF  (3  to  30  kHz)  and  LF  (30  to  300  kHz) 
bands  is  essentially  controlled  by  the  D  region  of  the  ionosphere.  Propagation  in 
the  ELF*  band  can  extend  into  the  E  and  F  regions.  Nuclear  weapon  effects  on  prop¬ 
agation  results  from  ionization  produced  by  the  detonation.  Usually  the  effects  are 
caused  by  free  electrons,  but  significant  absorption  may  result  from  ions  produced 
in  the  lower  part  of  the  D  region. 

Computer  codes  have  been  developed  by  the  Department  of  Commerce  (Office  of 
Telecommunications),  the  Department  of  the  Navy  (Naval  Ocean  Systems  Center),  and 
others  to  calculate  ELF,  VLF,  and  LF  propagation  for  both  natural  and  disturbed  condi¬ 
tions.  Generally,  these  codes  require  electron  and  ion  density  profiles  as  input  and 
perform  numerical  solutions  for  either  a  mode  or  ray-hop  sum  to  obtain  the  electric 
field  strength.  The  several  codes  are  characterized  by  whether  or  not  they  include 
the  geomagnetic  field,  the  method  used  to  include  earth's  curvature,  and  the  numerical 
methods  used.  An  important  assumption  often  used  in  the  calculational  models  is  that 
the  ionosphere  is  horizontally  stratified  (or  at  most  only  slowly  varying  between 
transmitter  and  receiver  terminals). 


* 

In  this  document  the  ELF  band  is  defined  as  3  to  300  Hz. 


The  WEDCOM  code  was  originally  developed  to  combine  nuclear  environments 
and  VLF  and  LF  propagation  models  in  a  single  computer  program.  The  first  version 
was  distributed  in  March  1969  and  provided  calculations  for  propagation  between 
vertically  polarized,  ground-based  antennas.  Input  requirements  (weapon  and  prop¬ 
agation  parameters)  were  those  generally  available  for  propagation  studies,  but 
were  made  sufficiently  flexible  so  that  the  code  could  be  used  in  nuclear  test  plan¬ 
ning,  data  analysis,  and  sensitivity  studies.  In  addition  to  outputs  required  to 
evaluate  received  field  strengths,  ionization  and  propagation  details  were  provided 
as  an  option. 

In  1970  a  second  version  of  the  code  (designated  WEDCOM  MB)  was  prepared 
for  use  in  calculating  effects  of  multiple  bursts.  In  addition,  ground  conductivity 
was  made  an  input  and  improvements  in  the  reflection  coefficient  calculations  were 
made.  In  1972  a  third  version  of  the  code  (designated  WEDCOM  III)  was  prepared  to 
include  calculations  of  VLF  horizontally  polarized  waves  and  to  include  VLF  height- 
gain  factors  for  elevated  transmitter  and  receiver  terminals.  In  addition,  improve¬ 
ments  in  the  ionization  models  and  code  structure  were  made. 

The  current  version  of  the  WEDCOM  code,  designated  WEDCOM  IV,  was  prepared 
in  1980  to  incorporate  improved  environment  models,  improved  VLF  and  LF  propagation 
models,  and  to  incorporate  a  model  for  ELF  propagation.  The  major  changes  to  environ¬ 
ment  model «  were: 

1.  Modification  of  D-region  chemistry  to  include  effect  of  cluster 
ions . 

2.  Modification  of  high-altitude  burst  phenomenology  affecting 
early-time  (minutes)  debris  location. 

3.  Addition  of  an  atmospheric  wind  model  affecting  late-time  (hours) 
debris  location. 

4.  Addition  of  new  prompt  ionization  sources  affecting  E-  and  lower 
F-region  ionization. 

The  VLF  and  LF  propagation  models  were  improved  as  follows: 

1.  Replacement  of  the  plane- slab  isotropic  ionosphere  with  an 
anisotropic  ionosphere  with  earth  curvature. 

2.  Provision  for  modeling  ground  conductivity  and  permittivity 
as  a  function  of  location. 

3.  Provision  for  elevated  and  arbitrarily  oriented  antennas  in 
the  LF  band . 
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4.  Addition  of  a  new  mode  search  model  in  the  VLF  band  to  insure 
obtaining  important  modes. 

5.  Addition  of  a  mode  conversion  program  in  the  VLF  band. 

1.2  LIMITATIONS  AND  DEFICIENCIES  OF  WEDCOM  IV 

The  major  limitations  and  deficiencies  of  WEDCOM  IV  identified  by  users 

include: 

1.  Link  evaluation  times  are  too  long.  While  the  propagation  models 
used  are  faster  running  than  the  detailed  models  developed  by 
NOSC,  they  still  require  too  much  evaluation  time  for  typical  user 
requirements . 

2.  Signal  variance  and  statistics  are  not  given. 

3.  Signal  processing  models  that  provide  either  bit  error  rates  or 
message  error  rates  are  not  given. 

4.  The  ambient  ionosphere  models  are  not  consistent  with  those 
used  by  NOSC  or  the  Office  of  Telecommunications  (which  also 
differ) . 

5.  Ambient  noise  models  are  not  provided  and  procedures  for 
determining  the  effect  of  nuclear  environments  on  received 
noise  are  not  included. 

6.  Signal- to- jamming  analysis  capability  is  not  included. 

7.  Code  structure  and  input  and  output  design  is  not  optimized 
for  system  analysis. 

As  noted  above  the  WEDCOM  IV  link  evaluation  times  are  significantly  longer 
than  for  previous  versions  of  the  code.  The  longer  running  times  are  due  to  includ¬ 
ing  anisotropy,  mode  search  and  mode  conversion  models,  and  an  improved  reflection 
coefficient  model.  The  longer  running  times  limit  use  of  the  code  in  studies  of  en¬ 
vironment,  propagation,  and  system  parameters.  Methods  for  reducing  link  evaluation 
times  are  discussed  in  Sections  2  through  6.  A  new  interpolation  procedure  for 
evaluating  VLF  and  LF  reflection  coefficients  is  described  that  reduces  computation 
time  while  still  providing  acceptable  accuracy.  A  new  VLF  mode  search  model  is 
described  that  is  faster  running  than  the  current  model  and  may  provide  significant 
improvements  in  running  time  for  undisturbed  or  slightly  disturbed  conditions. 

Other  model  improvements  that  can  be  included  in  WEDCOM  are  also  described  in  Sec¬ 
tions  2  through  4 . 
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A  new  ELF  propagation  model  developed  at  the  University  of  California,  San 
Diego  by  H.G.  Booker  that  requires  much  less  computation  time  than  the  current 
WEDCOM  model  is  described  in  Section  5. 

Some  reduction  in  computation  time  and  improvements  in  input  and  output 
formats  can  be  obtained  by  changes  in  code  structure.  Code  structure  considerations 
are  discussed  in  Section  6. 

The  deficiencies  described  in  items  2  through  4  have  been  studied  by  DNA 
contractors  and  models  are  available  that  can  be  included  in  a  WEDCOM  revision.  The 
deficiencies  in  items  5  through  7  are  largely  due  to  WEDCOM  being  a  single- link  code 
and  are  best  corrected  in  a  multilink,  multiburst  simulation  code  (eg,  SIMBAL) . 


► .  - 

*y- 
f  *'♦ 
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SECTION  2 

LF  MODEL  IMPROVEMENTS 


2.1  INTRODUCTION 

The  WEDCOM  IV  LF  propagation  model  was  revised  for  use  in  studies  for  the 
Naval  Electronic  Systems  Command  (Reference  2) .  The  revision  was  necessitated  by 
task  requirements  involving  high  (MF)  system  frequencies,  and  balloon  borne  antennas. 
The  model  improvements  were  concerned  with: 

•  A  more  efficient  and  simpler  program  structure. 

•  Better  treatment  of  tropospheric  refraction. 

•  Improved  convergence  factor. 

•  Improved  skywave  definition,  especially  for  elevated  antennas. 

•  Elimination  of  redundant  calculations. 

•  More  accurate  ionospheric  reflection  coefficient  computations. 

•  Improved  diffraction  models. 

•  Better  estimates  of  the  reflection  altitude. 

Elimination  of  redundant  computations  in  the  WEDCOM  IV  model  resulted  in  reduced 
computation  time.  However,  the  other  model  improvements  required  increased  computa¬ 
tion  time.  Thus,  a  computation  time  comparison  of  the  revised  model  with  WEDCOM  IV 
would  have  mixed  results.  For  example,  low-altitude  antenna,  low-frequency  cases 
with  a  uniform  ionosphere  over  the  propagation  path  would  execute  much  faster  with 
the  revised  model  (ie,  compared  to  WEDCOM  IV) ,  but  increased  computation  times  would 
probably  be  required  for  high-altitude  antenna,  high-frequency  cases  with  a  nonuniform 
ionosphere  over  the  path.  For  the  purposes  of  normal  WEDCOM  usage  the  high  accuracy 
requirements,  the  MF  cases,  and  high  balloon  altitude  antennas  are  not  important. 

Thus,  a  subset  of  the  above  improvements  can  be  incorporated  into  a  revised  WEDCOM 
model  that  will  increase  accuracy  while  not  requiring  a  significant  increase  in 
computational  time. 

Additional  improvements  to  the  LF  models  have  been  developed  for  the  present 
effort  that  will  further  improve  the  prediction  accuracy  and  computational  efficiency. 
Also,  several  new  and  revised  models  were  prepared  for  the  SIMBAL  code  (Reference  3) 
developed  for  DNA  and  DCA  that  can  be  used  in  the  WEDCOM  code  with  minor  revisirns. 

The  following  model  improvements  indicate  their  source: 


•  Analytic  reflection  coefficient  algorithm  (new,  S1MBAL) . 

•  Revised  LF  model  formulation  which  segments  the  propagation  path 
such  that  fewer  ionospheric  reflection  coefficient  computations 
are  required  and  improved  logic  implemented  to  eliminate  repeti¬ 
tive  computations  (new). 

•  An  improved  world  surface  conductivity  model  (SIMBAL) . 

•  An  improved  groundwave  computation  (SIMBAL) . 

•  A  more  efficient  (and  accurate)  diffraction  computation 
(new,  SIMBAL). 

A  major  cause  of  the  long  WEDCOM  IV  computation  times  is  the  ionospheric 
reflection  coefficient  computation.  Section  3  outlines  a  new  analytic  procedure  to 
reduce  this  requirement.  A  summary  of  the  remaining  LF  model  improvements  listed 
above  is  described  below.  Section  2.2  describes  the  new  procedure  for  defining  cal¬ 
culation  increments  along  the  propagation  path.  Section  2.3  describes  the  revised 
world  surface  conductivity  model.  Section  2.4  describes  the  revised  skywave  defini¬ 
tion.  An  improved  groundwave  computation  is  described  in  Section  2.5.  Sections  2.6 
through  2.8  briefly  describe  the  revised  LF  models  for  diffraction,  tropospheric  re¬ 
fraction,  and  ionospheric  convergence. 

2.2  LF  MODEL  FORMULATION 

The  WEDCOM  LF  model  presently  computes  ionospheric  electron- ion  profiles 
separated  by  equal  segments  along  the  propagation  path.  A  ground  conductivity  is 
specified  for  each  segment.  Adjacent  ionospheric  profiles  are  linearly  interpolated 
at  each  skywave  ionospheric  reflection  point  followed  by  the  reflection  coefficient 
computation.  Thus,  the  number  of  reflection  coefficient  computations  can  be  large, 
and  very  repetitive  for  uniform  profiles  where  anisotropic  effects  have  not  changed 
much.  A  second  area  of  large  computation  time  usage  is  the  number  of  repetitive 
diffraction  calculations  that  are  made  where  the  ground  constants  and  ionospheric 
reflection  characteristics  have  not  changed  from  one  skywave  hop  to  the  next  or  where 
they  were  computed  for  an  earlier  skywave  hop. 

Figure  1  illustrates  the  method  of  ground  and  ionospheric  path  segmentation. 
Ground  segment  endpoints  are  chosen  so  that  the  segments  are  of  equal  length  and 
segment  midpoints  occur  at  the  transmitter  and  receiver  locations.  Ground  conduc¬ 
tivity  calculations  are  made  at  the  ground  segment  midpoints.  Ionospheric  profiles 
(vertical  profile  of  ionization  and  collision  frequencies)  are  determined  at  loca¬ 
tions  chosen  as  described  in  Section  6.  The  spacing  between  ionospheric  profiles 
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Figure  1.  Illustration  of  LF  ground  and  ionospheric 
path  segmentation. 


can  differ  and  the  ionospheric  profile  locations  are  independent  of  the  ground  seg¬ 
ment  locations .  Ionospheric  segments  are  chosen  so  that  ionospheric  segment  end¬ 
points  occur  at  the  transmitter  and  receiver  and  at  the  ground  segment  endpoint 
nearest  the  midpoint  between  ionospheric  profile  points. 

Reflection  coefficient  computations  are  made  at  the  first  ground  segment 
midpoint  within  an  ionospheric  segment,  and  thereafter  where  the  anisotropic  ef¬ 
fects  change  significantly  at  a  new  ground  segment  midpoint.  This  criteria  continues 
until  the  next  ionospheric  segment  has  been  reached.  The  end  result  is  a  number  (NR(,) 
of  reflection  coefficient  segments.  The  number  of  reflection  coefficient  segments 
can  be  large  without  a  noticeably  increased  computation  time,  but  will  provide: 

•  Realistic  ground  conductivity  changes  along  the  propagation  path. 

•  Convenient  ionosphere  segmentation  where  each  segment  represents 
a  precomputed  ionospheric  profile. 

•  Anisotropic  ionosphere  sensitivities  (tested  at  each  ground 
segment  midpoint)  within  an  ionospheric  segment . 

Figure  1  illustrates  this  computation  criterion  where  no  additional  computations 
occur  due  to  ionospheric  anisotropy  (N^  is  equal  to  4) . 

Reflection  altitudes  at  the  locations  used  to  compute  ionospheric  profiles 
are  computed  as  described  in  Section  3.  Then  the  skywave  geometry  is  computed  by 
linearly  interpolating  the  reflection  altitudes.  This  procedure  is  not  conceptually 
different  than  in  WEDCOM  IV  or  the  LF  model  revised  for  the  balloon  studies.  Now, 
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however,  the  number  of  ionospheric  reflections  and  their  incidence  angles  are  kept 
track  of  for  all  skywaves  within  each  reflection  coefficient  segment.  The  steepest 
and  most  oblique  angle  of  incidence  within  any  one  of  the  segments  is  known.  Thus, 
the  number  of  reflection  coefficient  computations  within  a  segment  can  be  minimized 
by  the  analytical  reflection  coefficient  representation  described  in  Section  3. 
Normally  within  a  reflection  coefficient,  one  or  two  (between  extreme  oblique  and 
steep  incidence)  computations  need  be  performed  for  most  cases.  Possibly  for  some 
selected  cases  a  few  more  computations  are  necessary  (for  example,  at  the  higher  LF 
frequencies  where  smaller  incidence  angle  spacing  is  required  for  accurate  interpola 
tion) .  Nevertheless,  the  number  of  computations  are  significantly  reduced  from  the 
WEDCOM  IV  version. 

The  revised  LF  model  developed  for  the  balloon  studies  addressed  the 
repetitive  nature  of  the  ionospheric  reflection  coefficient  and  diffraction  calcula¬ 
tions.  Logical  tests  were  added  that  used  the  computations  from  a  prior  skywave-hop 
if  the  present  skywave-hop  has  the  same  geometrical  and/or  ground  characteristics. 
Section  3  describes  an  improved  criteria  for  testing  anisotropy  sensitivity.  The 
present  effort  identified  additional  internal  repetitive  computations  within  the  dif 
fraction  model  which  were  eliminated  through  logical  tests.  The  diffraction  model 
improvements  are  discussed  in  Section  2.5. 

2.3  WORLD  SURFACE  CONDUCTIVITY 

The  ground  conductivity,  a,  and  relative  permittivity,  er>  strongly  affect 
the  amount  of  reflected  electromagnetic  wave  losses  from  the  earth's  surface.  A 
large  range  of  values  exist — typical  values  are: 


a 

er 

Conductivity 

Relative 

Location 

mho/m 

Dielectric  Constant 

Sea  Water 

4 

81 

Fresh  Water 

nr3 

81 

Wet  Earth 

ID’3 

18 

Dry  Earth 

ID'5 

5 

The  importance  of  the  seasonal  water  content  of  the  soil  can  be  noted. 

Fortunately,  a  detailed  conductivity  map  is  not  required  for  ELF/VLF/LF 
propagation.  Indeed,  an  exact  map  may  lead  to  inaccuracies  as  small  regions  of 
conductivity,  differing  from  conductivities  of  larger  surrounding  regions,  should 
not  significantly  affect  propagation  and  should  not  be  modeled.  Thus,  useful 
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conductivity  maps  should  be  gross,  with  only  a  few— less  than  10—  (cr,er)  pairs 
representing  the  world  conductivity. 

The  WEDCOM  IV  ground  conductivity  model  (Subroutine  CONMAP)  is  a  modified 
version  of  the  Institute  of  Telecommunication  Sciences  world  conductivity  map  (Ref¬ 
erence  4).  Here  a  Fourier  analysis/polynomial  interpolation  procedure  was  applied 
to  a  world  conductivity  map.  The  model  uses  the  Fourier  coefficient/polynomial 
representation  to  determine  a  parameter  (here  called  y)  at  a  given  location.  The 
conductivity  (a)  and  relative  dielectric  constant  (e^)  values  are  determined  from 
the  value  of  y  (as  seen  in  Figure  2  for  a) .  Note  the  continuous  variation  of  a 
values.  This  part  of  the  model  is  inconsistent  with  the  "gross"  requirement  for  low- 
frequency  propagation,  and  is  especially  inaccurate  near  water- land  transitions. 
Thus,  not  only  is  an  inaccuracy  introduced,  but  the  increased  number  of  possible 
(0, er)  values  limits  the  effect  of  more  efficient  (faster  computational)  propagation 
modeling. 

It  is  proposed  to  replace  routine  CONMAP  with  the  model  used  in  the  SIMBAL 
code  (routine  RDCS,  Reference  3).  The  SIMBAL  model  utilizes  a  worldwide  tabulated 
data  and  provides  the  gross  (cr,er)  map  desired.  It  is  a  simple  "table  lookup"  with 
no  computational  elegance  as  done  in  CONMAP.  The  outputs  are  always  discrete— not 
continuous.  There  are  only  six  (a,er)  pairs  output.  The  only  disadvantage  in  using 
the  SIMBAL  model  is  the  data  storage  (5.2  K  versus  the  1.5  K  in  CONMAP).  However, 
this  should  be  outweighed  by  the  improved  predictions.  Although  smaller  longitudinal 
increments  are  possible  the  minimum  latitudinal  variation  is  5  degrees  or  approxi¬ 
mately  500  km.  Thus  it  is  difficult  to  justify  ground  segments  less  than  approx¬ 
imately  250  to  300  km. 

2.4  SKYWAVE 

Figure  3  illustrates  the  ray-path  geometry  types  found  in  the  geometric- 
optic  model  revised  for  the  balloon  study.  A  description  of  the  WEDCOM  IV  geometry 
will  not  be  repeated  here  except  to  note  that  it  is  a  subset  of  these  geometries. 
While  only  a  single  ionospheric  reflection  is  shown,  extension  to  multiple  hops  is 
straightforward.  Figure  3a  illustrates  four  "Earth-Detached"  ray  paths  (designated 
as  "WG"  for  the  classical  "Whispering  Gallery"  phenomena) .  They  are  so  termed  as 
they  never  touch  the  earth's  surface.  Because  it  was  awkward  to  draw,  ray-path  "2" 
is  shown  to  coincide  with  ray-path  "1."  Actually,  ray-path  "1"  is  tangent  to  the 
earth's  surface  near  each  end,  and  these  tangent  points  lie  at  lower  altitudes  than 
the  antenna  altitudes.  A  fifth  ray-path  could  exist  that  (for  a  single  hop)  would 
also  not  touch  the  earth — except  by  the  path  extension  beyond  the  RCVR  or  XMTR 


Figure  2.  Conductivity  prediction  as  a  function  of  the  parameter  y  for 
Subroutine  CONMAP  used  in  WEDCOM  IV. 

locations.  This  path  is  shoun  in  Figure  3b  (ray- path  "5")  and  is  not  a  WG  ray-path. 
All  five  types  do  not  coexist  for  any  given  number  of  ionospheric  reflections,  but  a 
combination  of  them  can.  Figure  3b  also  shows  the  remaining  (E-E)  ray-path  types. 
These  ray- paths  will  reflect  off  the  earth’s  surface.  The  elevated  antenna -to -ground 
antenna  case  (E-G)  (Figure  3c)  illustrates  an  earth-surface  diffracted  path  and  a 
single  earth  reflected  ray.  For  the  diffracted  case  the  ray-path  is  tangent  at  the 
earth's  surface  where  it  diffracts  the  remaining  range  to  the  RCVR.  The  third 
geometry  mode,  ground  antenna- to -ground  antenna  propagation  paths,  is  illustrated 
in  Figure  3d. 

Figure  3b  illustrates  the  5,  6,  7,  and  8  ray-paths  types.  Note  in  Figures 
3c  and  3d  how  these  ray-path  types  coincide  as  one  or  both  antennas  are  on  the  ground 
Table  1  lists  the  optic-region  foreground  factors  for  all  eight  ray-path  cases.  For 
diffracted  cases  (here  the  WG  ray-paths  are  not  applicable)  the  foreground  factors 
all  correspond  to  the  1+R  geometric-optic  case,  as  the  ray-path  types  6,  7,  and  8 
coincide  with  ray-path  type  "5." 
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XMTR 

RCVR 


(a)  BOTH  ANTENNAS  ELEVATED  (WG) 

(Whispering  Gallery  or  Earth-Detached) 


(b)  BOTH  ANTENNAS  ELEVATED  (E-E) 


ONE  ELEVATED  AND  ONE  GROUND  ANTENNA  (E-G,  G-E) 


5,6, 7,8 


(d)  BOTH  GROUND  ANTENNAS  (G-G) 

Illustrating  the  possible  single  Ionospheric  reflecti 
ray- path  geometry. 


Table  1.  Foreground  factors  for  geometric  optic 
ray-paths  (nondiffracted). 


Antenna 

Case 


Ray-Path  Type 


•  R  is  defined  as  the  phase-corrected  ground 
reflection  coefficient. 

•  Arrows  show  the  coincidence  of  ray-path  types. 


Figure  4  illustrates  typical  daytime  propagation  ranges  versus  transmitter 

antenna  elevation  angle  for  ray-paths  possessing  2,3 - 7  ionospheric  reflections. 

The  ray-path  types  defined  in  Figure  3  are  indicated.  The  "peak"  values  correspond 
to  where  the  ray-path  is  tangent  to  the  earth's  surface.  The  dotted  lines  denote 
approximate  regions  where  diffraction  modeling  must  be  used,  and  the  approximate 
Earth-Detached  (WG)  and  geometric-optics  regions.  In  the  diffraction  region,  the 
ray-path  types  are  grouped  (6,4),  (8,1),  (7,3),  and  (5,2).  The  nine-Mm  range  is 
selected  as  an  example  of  how  this  figure  could  be  used  to  predict  the  significant 
ray- paths  and  how  they  are  modeled.  For  this  c/se  there  are  no  WG  ray-paths  as  they 
are  within  the  diffraction  regions.  The  following  lists  the  significant  ray-paths 
(diffraction  modeled  ray-paths  in  parenthesis): 


Ionospheric 

Reflections 


Ray-Path 


5, 6, 7, 8 
5, 6, 7, 8 

(5.2) , (6,4), (7,3), (8,1) 

(5.2) *, (6,4)*, (7,3)*, (8,1) 


8 


The  asterik  (*)  ray-paths  diffract  around  the  earth's  surface  to  achieve  the  9-Mm 
range. 

The  ray-path  characteristics  are  found  in  a  manner  similar  to  that  usea  in 

WEDCOM  IV.  First,  the  ground  elevation  angle  is  set  to  zero  for  a  given  hop  number. 

The  range  is  computed.  If  this  range  is  shorter  than  necessary,  the  hop  number  is 
increased  until  the  range  is  equal  to  or  exceeds  the  required  range.  A  very  fast 
Newton  iteration  algorithm  is  then  used  to  find  the  ray-path  geometry.  At  least  one 

(ie,  for  the  G-G  mode)  earth  diffracted  ray-path  is  defined,  and  several  are  defined 

for  elevated  antenna  cases.  At  least  two  (G-G  mode)  additional  ray-paths  with  in¬ 
creased  ionospheric  reflections  are  also  defined.  Here  again,  for  elevated  antennas 
(ie,  for  the  E-G,  G-E,  and  E-E  modes)  this  results  in  a  number  of  individual  ray- 
paths. 

2.5  GROUNDWAVE 

The  WEDCOM  IV  groundwave  field  computation  described  in  Reference  1  is 
strictly  applicable  to  a  uniform  ground  conductivity  propagation  path.  For  this 
case  only  a  single  computation  using  that  procedure  is  required.  For  a  variable 
ground  conductivity,  WEDCOM  IV  took  an  average  value  of  path  ground  conductivity 
and  again  applied  the  same  procedure.  The  inadequacy  of  this  method  is  obvious  for 
paths  propagating  partially  over  sea  water  with  the  remainder  being  low-ground 
conductivity. 

A  combination  of  physical  reasoning  and  heuristic  argument  is  described  by 
Millington  (Reference  5)  whereby  the  computational  procedure  of  Reference  1  is  ap¬ 
plied  iteratively  over  the  variable  path  conductivity  segments.  The  Millington  pro¬ 
cedure  has  shown  amazing  agreement  with  more  sophisticated  computation  procedures 
(References  6  and  7).  For  example  the  large  field  strength  recovery  crossing  from 
land  to  sea  is  predicted  by  Millington's  procedure.  Moreover,  the  signal  phase 
prediction  is  very  accurate,  a  surprising  phenomena  as  Millington  made  no  considera¬ 
tion  of  the  phase  change  at  the  boundary  of  a  conductivity  change. 

Eckersley  (Reference  8)  is  first  credited  with  the  procedure  of  utilizing 
uniform  ground  conductivity  computations  to  piece  together  a  variable  conductivity 
path  prediction.  Figure  5  demonstrates  the  Eckersley  procedure  for  a  three- segment 
path.  The  segment  including  the  transmitter  has  a  ground  conductivity  pair  (o^.e^), 
the  center  segment  (a 2,^2)’  arK*  t^ie  segment  including  the  receiver  (a^,c^).  The 
three  (a,e)  pair  curves,  Sj,  S2,  and  S^,  are  calculated  by  the  method  of  Reference  1. 
The  fat  solid  line  in  Figure  5  represents  the  predicted  received  field  strength 
using  the  relationships 
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Eckersley's  method  does  not  obey  the  principle  of  reciprocity  (except  for  a  homo¬ 
geneous  path,  a  degenerate  case) . 

Millington's  method  does  satisfy  the  principle  of  reciprocity  as  the  pre¬ 
dicted  received  signal  is  the  geometric  mean  of  Eckersley's  method  in  Figure  3a 
3.  b 

(S^)  and  Figure  3b  (S^)  where  the  transmitter  and  receiver  locations  have  been 
interchanged.  That  is,  the  predicted  received  field  strength  by  Millington's  pro¬ 
cedure  is 
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was  defined  above  and 
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This  example  had  three  segments.  Millington  generalized  the  procedure  to  include 
N  segments: 
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The  parenthesis  groupings  in  Equation  6  correspond  to  (a.e^)  pairs.  The  implementa¬ 
tion  of  Equation  6  will  incorporate  the  following  items  for  computational  efficiency 
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•  A  single  computation  will  be  used  if  the  path  is  homogeneous. 

•  Where  adjacent  (o,er)  pairs  are  nearly  identical  they  will  be 
merged  together  so  that  only  the  minimum  parenthetical  groupings 
will  be  evaluated. 

•  The  parenthetical  computations  are  added  together  (rather  than 
multiplied)  as  the  £tt(S^)  is  computed. 

2.6  EARTH  SURFACE  DIFFRACTION,  REFLECTION 

Section  2.4  described  the  skywave  geometry  for  the  regions  where  geometric- 
optics  is  valid,  where  diffraction  is  the  dominant  characteristic,  or  in  the  inter¬ 
mediate  (caustic)  region  (see  Figure  6).  The  model  for  analyzing  the  effects  due  to 
the  earth's  surface  is  now  examined.  The  detailed  expressions  used  by  WEDCOM  IV  for 
the  ground  Fresnel  reflection  coefficients  and  the  diffraction  are  found  in  Reference 
1  and  won't  be  repeated  here. 

The  caustic  region  (see  Figure  6b)  has  been  a  source  of  model  inconsis¬ 
tency.  WEDCOM  IV  switched  from  the  geometric  optics  model  to  the  diffraction  model 
when  the  ground  elevation  angle,  e,  became  less  than  0.03  radians.  This  switch 
point  did  not  provide  accurate  predictions  for  many  cases.  An  improved  switch 
criteria  is  addressed  here. 

The  WEDCOM  IV  diffraction  model  is  a  major  user  of  computation  time,  sec¬ 
ond  only  to  the  ionospheric  reflection  coefficient  model.  There  are  many  repetitive 
computations  in  the  WEDCOM  IV  model,  especially  when  the  ground  conductivity  is 
unchanged  for  large  parts  of  the  propagation  path,  or  when  path  segments  had  the 
same  conductivity,  but  were  separated  by  segments  of  differing  conductivity.  This 
repetition  was  eliminated  in  the  newer  diffraction  models  thereby  reducing  (for  many 
cases)  the  computation  times  by  several  factors. 

Table  2  lists  the  cases  for  which  both  the  geometric-optic  and  diffrac¬ 
tion  models  were  exercised.  The  caustic  region  is  illustrated  in  Figures  7  through 
10  for  an  ionosphere-to- ionosphere  case,  and  Figures  11  through  14  for  ionosphere- 
to-antenna  cases.  A  smooth  transition  from  the  diffraction  to  geometric  optics 
models  is  not  always  possible.  This  is  often  true  for  the  ionosphere-to-ionosphere 
cases  and  is  shown  in  Figures  7  through  10.  The  new  model  logic  will  select  the 
lowest  magnitude  value  between  the  two  models  which  minimizes  the  abrupt  transition. 
The  abrupt  transition  is  due  in  part  to  the  diffraction  model  effectively  computing 
a  "1+R"  value  while  the  geometric-optic  model  assumes  only  the  Fresnel  reflection 
coefficient  "R"  ray.  Improved  transition  computation  involves  a  refinement  of  the 


23 


XMTR , ION 


RCVR , ION 


EARTH'S  SURFACE 


(a)  Geometrical  Optics  Region 


XMTR, ION 


XMTR-TRANSMITTER  ANTENNA  LOCATION 
RCVR-RECEIVER  ANTENNA  LOCATION 
ION- IONOSPHERIC  REFLECTION  LOCATION 


RCVR, ION 


EARTH'S  SURFACE 


(b)  Caustic  Region 


XMTR, ION 


EARTH'S  SURFACE 


(c)  Diffraction  Region 


RCVR, ION 


Figure  6.  Illustrating  how  the  geometric-optic  region  ray 

paths,  a  and  b,  converge  in  the  diffraction  region, 


FOREGROUND 


Table  2.  Matrix  of  cases  exercised  by  the  geometric-optics  and 
diffraction  models. 


Ground 

Reference 

Reference 

Frequency 

(o»e) 

Altitude  1 

Altitude  2 

(mhos/m) 

(km) 

_ (km) _ 

GROUND  ELEVATIO-  -  ■  , radians) 

Figure  7.  Comparison  of  geometric-optic  and  diffraction  models  of 
foreground  factor  in  caustic  region  for  ionosphere-to- 
ionosphere  (Hr  =80  km),  sea  water,  and  vertical 
polarization. 
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FOREGROUND  FACTOR  MAGNITUDE 


-0.03  -0.018  -0.009  0  0.009  0.018  0.027  0.036 


GROUND  ELEVATION  ANGLE  (radians) 


Figure  9.  Comparison  of  geometric-optic  and  diffraction  models  of  foreground 
factor  in  caustic  region  for  ionosphere- to- ionosphere  (Hr =80  km), 
land,  and  horizontal  polarization. 
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GEOMETRIC-OPTIC 

MODEL 


GROUND  ELEVATION  ANGLE  (radians) 

Figure  10.  Comparison  of  geometric -optic  and  diffraction  models  of  foreground 
factor  in  caustic  region  for  ionosphere-to-ionosphere  (Hr  =  80  km), 
land,  and  horizontal  polarization. 
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Figure  11.  Comparison  of  geometric-optics  and  diffraction  models  of  foreground 
factor  in  caustic  region  for  ionosphere-to-antenna  case  (Hr  =  80  km, 
antenna  altitude  =0.0  km),  land,  and  vertical  polarization. 
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Figure  13.  Comparison  of  geometric  optics  end  diffraction  models  of  foreground 
factor  In  caustic  region  for  Ionosphere- to-antenna  case  (Hr  -  80  km, 
antenna  altitude- 10  km),  land,  and  vertical  polarization. 
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skywave  definitions  (Section  2.4)  which  will  be  addressed  when  this  model  is  imple¬ 
mented  into  WEDCOM.  For  the  most  part,  however,  the  abruptness  is  due  to  the  sep¬ 
arate  consideration  of  the  direct  "1,"  and  reflected  "R"  rays  in  the  geometric-optic 
region. 

The  transition  between  the  two  models  is  seen  to  be  much  smoother  for  the 
ionosphere- to- antenna  cases  (Figures  11  through  14).  As  mentioned  above,  the  WEDCOM 
IV  model  presently  switches  at  1.72  degrees  (0.03  radians).  An  improved  value  has 
been  selected  based  on  the  Table  2  matrix  of  computations.  This  value  is  5  degrees 
which  will  have  a  tendency  to  increase  the  signal  strength  from  prior  model  versions. 

2.7  TROPOSPHERIC  REFRACTION 

The  troposphere  has  been  long  recognized  as  affecting  radio  waves  launched 
at  low  elevation  angles  at  VHF  frequencies  and  above.  Groundwave  signals  at  LF 
frequencies  (for  example,  Loran)  are  also  affected  by  the  troposphere  refraction  ef¬ 
fects.  Agreement  between  signal  measurements  and  relatively  simple  analytic  predic¬ 
tion  techniques  (for  example,  equivalent  earth  radius)  have  made  this  propagation 
phenomena  well  understood.  Less  clear  are  the  effects  of  the  troposphere  on  prop¬ 
agation  when  energy  escapes  the  troposphere  and  is  reflected  by  the  ionospheric 
region.  It  is  the  latter  case  that  is  of  interest  in  LF  skywave  propagation. 

WEDCOM  IV  utilized  an  equivalent  earth's  radius  in  the  LF  model  to  account 
for  tropospheric  refraction.  The  VLF  model  had  no  refraction  modeling  and  the  com¬ 
parison  of  LF  and  VLF  model  predictions  were  not  favorable  due  in  part  to  this  reason. 
Tropospheric  refraction  at  VLF  and  LF  was  addressed  in  Reference  2.  A  summary  of 
results  is  presented  here. 

The  troposphere  had  almost  no  effect  at  VLF,  and  little  effect  for  the 
lower  LF  frequencies,  especially  when  ionospheric  reflection  was  strong.  There  was 
a  noticeable  effect  for  the  higher  LF  cases,  especially  when  the  ionospheric  reflec¬ 
tion  was  weak  where  the  differential  reradiated  field  in  the  troposphere  was  compa¬ 
rable  to  the  reradiated  field  at  the  ionospheric  reflection  altitudes.  Even  for  the 
cases  where  tropospheric  refraction  is  important,  however,  the  equivalent  earth's 
radius  model  (ie,  4/3  a,  where  a  is  the  actual  earth's  radius)  grossly  overestimated 
the  effect  on  skywave  propagation. 

The  convenience  of  a  ray  theory  treatment  (such  as  the  equivalent  earth's 
radius  technique)  of  tropospheric  refraction  is  obvious.  An  accurate  ray-trace 
model  is  inconsistent  with  the  objectives  of  WEDCOM.  The  new  LF  tropospheric  refrac¬ 
tion  model  retains  an  equivalent  earth's  radius,  but  with  different  values  empirically 


selected  in  Reference  2.  The  model  characteristics  are  as  follows: 


1.  An  equivalent  earth  radius  of  4/3  a  will  be  used  for  all  dif¬ 
fracted  rays  and  groundwaves  for  frequencies  near  65  kHz  and 
above.  The  equivalent  earth's  radius  will  decrease  linearly 
with  frequency  from  ag  =  8485  km  at  65  kHz  to  ag  =  6364  km  at 
45  kHz. 

2.  The  earth  detached  modes  are  not  affected  by  the  troposphere 
and  will  use  normal  earth  curvature  (ie,  a  =  6364  km) . 

3.  Skywaves  that  reflect  from  the  earth's  surface  will  use 

ag  =  7000  km  for  65  kHz  and  above.  A  linear  decrease  (with 
frequency)  to  ag  =  6364  km  at  45  kHz  will  be  made. 

4.  Actual  earth's  radius  will  be  used  for  all  computations  below 
45  kHz. 

2.8  CONVERGENCE  FACTOR 

The  ionospheric  convergence  coefficient  is  a  geometric  factor  that  ac¬ 
counts  for  the  convergence  of  rays  reflected  from  an  idealized  spherical  reflector. 

Two  important  approximations  are  made  in  the  definition  of  the  convergence  coef¬ 
ficient.  One  is  that  the  geometric  optics  assumption  applies,  and  the  other  is  that 
energy  is  reflected  from  a  thin,  well-defined  region  in  the  ionosphere.  However, 
the  reflection  region  is  not  always  well  defined  (Reference  9  addressed  this  case) . 
Also,  when  the  ray  approaches  the  caustic  point  the  geometrically  defined  conver¬ 
gence  factor  becomes  infinite  since  in  the  vicinity  of  the  caustic  the  geometric 
optics  assumption  fails. 

Wait  (Reference  10)  has  developed  a  correction  factor  based  on  wave  theory 
which  cancels  the  infinity  and,  in  general,  limits  the  convergence  coefficient  in  the 
vicinity  of  the  caustic  to  about  a  factor  of  two.  Reference  10  also  contains  an  ap¬ 
proximation  to  extend  the  definition  of  the  convergence  coefficient  beyond  the  caustic 
point  for  surface  antennas . 

When  the  antennas  are  elevated,  the  geometric  convergence  factor  can  be 
very  large  and  remains  large  even  after  applying  the  wave  theory  correction. 
Straightforward  application  of  the  geometric  definition  plus  wave  corrections  leads 
to  high  peak  field  values  at  specific  locations  that  depend  critically  on  path 
distance  and  ionospheric  reflection  altitude.  The  sensitivity  of  the  convergence 
coefficient  to  reflection  altitude  was  addressed  in  Reference  9. 
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Figure  15  shows  the  corrected  convergence  coefficient  as  a  function  of  path 
distance  for  transmitter  and  receiver  altitudes  of  25  km  and  for  parametric  reflec¬ 
tion  altitudes.  Note  that  there  is  little  sensitivity  to  reflection  altitude  at  dis¬ 
tances  less  than  the  distance  to  the  caustic,  but  significant  variation  at  greater 
distances.  In  particular,  the  distance  to  produce  the  peak  in  the  convergence  factor 
varies  with  the  reflection  altitude.  Since  the  reflection  region  may  be  several 
kilometers  thick,  a  weighted  average  over  the  reflection  region  is  used  to  estimate 
the  convergence  coefficient. 

The  new  LF  convergence  coefficient  model  in  the  region  beyond  but  near  the 
caustic  limits  the  geometric  convergence  coefficient  to  be  no  greater  than  three. 

This  removes  the  large  peaks  in  the  field  strength  that  occur  for  particular  distance- 
reflection  altitude  combinations,  and  causes  the  total  field  to  vary  relatively 
smoothly  with  distance  as  the  receiver  moves  away  from  the  caustic  point  and  the 
geometric-optics  assumption  again  applies. 
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Figure  15.  Convergence  coefficient  sensitivity  to  reflection  altitude. 
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SECTION  3 


IONOSPHERIC  REFLECTION  ALTITUDE 
AND  REFLECTION  COEFFICIENT  MODELS 


3.1  INTRODUCTION 

A  reflection  altitude  calculation  is  required  in  both  the  LF  and  VLF  prop¬ 
agation  models.  A  description  of  an  improved  procedure  for  determining  a  reflection 
altitude  algorithm  is  described  in  Section  3.2.  Reflection  coefficient  calculations 
are  made  at  ground  segment  midpoints  as  described  in  Section  2.2.  Section  3.3  de¬ 
scribes  improved  criteria  for  determining  when  anisotropic  effects  change  sufficiently 
to  require  reflection  coefficient  calculations. 


The  ionospheric  reflection  coefficient  computation  is  the  single  largest 
user  of  computer  time  of  all  the  WEDCOM  models.  Both  the  waveguide  and  geometric- 
optic  solutions  require  this  computation  to  be  repeated  numerous  times.  It  is  thus 
desirable  to  eliminate  unnecessary  computations  or  scale  existing  computations  when¬ 
ever  possible.  Analytic  forms  of  the  reflection  coefficient  are  described  here  that 
will  provide  the  proper  balance  between  computing  time  and  accuracy.  Before  pro¬ 
ceeding,  however,  it  is  appropriate  to  review  the  background  and  nature  of  the  ionos¬ 
pheric  reflection  coefficient  computation. 


The  ionospheric  reflection  coefficients  can  be  represented  as  a  four-com¬ 
ponent  matrix,  ie, 


where 


R  = 


IIR||  iRll 


.«Ri  iRX 


(7) 


R  =  reflection  coefficient  matrix 

ii  =  polarization  parallel  to  the  plane  of  incidence 
l  =  polarization  perpendicular  to  the  plane  of  incidence. 

The  first  subscript  applies  to  the  polarization  of  the  incident  wave,  and  the  second 
subscript  applies  to  the  polarization  of  the  reflected  wave. 


The  coupling  coefficients,  AR>.  and  i,Ra,  may  be  as  large  or  larger  than  the 
primary  coefficients  „R,,  and  AR A  for  normal  nighttime,  but  they  are  typically  much 
smaller  than  the  primary  coefficients  for  daytime  conditions.  For  a  normal  ionos¬ 
phere,  the  coefficients  are  a  function  of  propagation  azimuth  and  the  earth's  magnetic 
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field.  A  nuclear  disturbance  results  in  lowering  the  reflection  region.  As  a  con¬ 
sequence  of  the  increased  electron-neutral  collision  frequency  in  the  reflecting 
region,  a  strong  disturbance  can  cause  the  coupling  coefficients  to  be  essentially 
zero. 

LF  and  VLF  waves  are  reflected  from  the  D-region.  Most  of  the  reflected 
energy  is  returned  from  the  60-  to  70-km  altitude  level  for  normal  daytime  ionos¬ 
pheres  and  from  80-  to  90-km  altitude  level  for  normal  nighttime  ionospheres.  In 
these  altitude  ranges,  particularly  for  nighttime  conditions,  the  electron-neutral 
collision  frequency  is  comparable  to  the  electron  gyro frequency,  and  the  earth's 
magnetic  field  has  a  significant  effect  on  the  reflection  coefficient. 

Presently,  the  WEDCOM  IV  code  (Reference  1)  employs  two  calculation  pro¬ 
cedures  for  determining  the  ionospheric  reflection  coefficients.  The  first  method 
is  an  iterative  integration  of  a  first-order  differential  equation.  It  is  applic¬ 
able  to  both  isotropic  and  anisotropic  ionospheres.  The  second  procedure  is  an  ap¬ 
plication  of  a  simple  recursive  technique  and  is  used  only  with  isotropic  ionospheres. 

The  iterative  integration  is  a  very  complex  calculation  and  its  efficiency 
depends  critically  on  well-defined  starting  (highest)  altitudes  and  stopping  (lowest) 
altitudes,  numerical  error  criteria,  and  integration  stepsize.  In  the  past,  this 
computation  was  economized  by  adjusting  these  constraints.  While  this  technique  did 
reduce  the  running  time,  it  remained  substantial.  Moreover,  for  some  applications 
the  "adjustments"  resulted  in  unacceptable  integration  accuracies.*  Setting  the 
integration  criteria  to  values  that  assure  acceptable  accuracy  for  all  cases  requires 
a  more  stringent  criteria  than  now  found  in  WEDCOM  IV. 

A  significant  reduction  in  computation  time  can  be  obtained  if  an  analytic 
form  for  the  reflection  coefficient  can  be  used  where  the  coefficients  are  determined 
from  only  a  few  reflection  coefficient  calculations.  This  has  been  done  in  WEDCOM 
for  isotropic  conditions.  Previous  attempts  to  find  useful  analytic  forms  for  aniso¬ 
tropic  conditions  have  only  been  partially  successful  in  that  they  resulted  in  poor 
modeling  over  part  of  the  region  of  interest.  Recent  work  described  in  Section  3.4 
has  lead  to  an  applicable  analytic  form  that  can  be  used  in  WEDCOM  for  both  isotropic 
and  anisotropic  conditions. 

_ 

Inaccuracies  tend  to  occur  as  the  frequency  is  increased,  for  the  more  steep 
incident  angles,  and  as  the  ionospheric  profile  becomes  lossy  over  large  regions. 
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3.2  REFLECTION  ALTITUDE 

WEDCOM  IV  presently  uses  a  simple  analytic  formulation  to  determine  the 
reflection  altitude,  HR.  It  is  a  simplified  derivative  of  the  more  accurate 
expression 


Re (cos' 


v  - 


2-0|x|/( 


,o  2  „  1/2 

(8  +  cos  £)  -  cos  £ 


(8) 


where  Re  means  "real  part  of,"  x  is  the  isotropic  susceptibility  (defined  in  the 
next  section)  evaluated  at  HR 


£  = 


tan'1  Im(x)/Re(x)J 


(9) 


and  0^  is  the  angle  of  incidence  at  the  ionosphere  (see  Reference  1) .  The  relation¬ 
ship  in  Equation  8  is  satisfied  at  a  single  altitude  (ie,  HR)  and,  thus  the  value 
of  HR  is  insensitive  to  signal  propagation  effects  below  or  above  this  altitude. 
Actually,  reflection  occurs  at  every  point  in  the  medium  and  to  compute  the  effective 
reflection  altitude  accurately,  the  contribution  for  each  altitude  needs  to  be  summed 
for  the  effects  on  the  signal  seen  below  the  reflection  region. 

Detailed  studies  were  made  in  References  2  and  3  which  indicated  that  HR, 
while  always  high  as  derived  from  Equation  8,  was  reasonably  accurate  when  the 
ionospheric  reflection  region  was  very  sharp.  HR  was  found  to  be  much  too  high  for 
some  disturbed  profiles  where  the  reflection  region  was  large  or  when  the  reflection 
process  was  weak.  This  is  illustrated  in  Figures  16  and  17  where  the  up  and  down¬ 
going  H-field  and  Poynting  vector  components  are  shown.  Note  the  sharp  reflection 
gradient  for  the  normal  nighttime.  For  the  studies  in  Reference  2  an  effective 

reflection  altitude,  HP,  was  selected  at  the  point  where  the  downgoing  H  component 

c  y 

(HYD)  is  one-half  the  total  reflected  value.  Thus,  while  HR  is  high  it  is  not  too 
excessive  for  normal  night.  For  the  disturbed  environment,  however,  Hg  is  consider¬ 
ably  below  Hr.  Figure  18  shows  the  value  of  HR  and  the  altitude  range  from  which  an 
effective  reflection  altitude  could  be  selected. 


The  determination  of  Hg  as  performed  above  is  not  computationally  practical 
for  the  purposes  of  WEDCOM.  An  analytic  form  such  as  Equation  7  is  still  required. 
The  value  of  HR  is  high  as  it  doesn't  actually  account  for  the  signal  absorption  at 
lower  altitudes.  Also  shown  on  Figure  18  is  a  corrected  value  of  11R,  HR(,,  which  was 
made  by  reducing  HR  by  H  ,  ie 


H, 


RC 


Hd  -  H 
R  X 


(10) 
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Figure  16.  H-field  component  magnitudes  and  Poynting  vector  components  in  vertical  direction  (z)  for 

normal  nighttime  conditions,  for  40  kHz  and  ei  =  80°.  x  is  the  direction  of  propagation  and 
y  is  the  direction  orthogonal  to  x,z  and  perpendicular  to  the  direction  of  propagation. 
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(a)  H-field  0>)  Poynting  vector 

Figure  17.  H-field  component  magnitudes  and  Poynting  vector  components  in  vertical  direction  (z)  for 

strong  gamma  nuclear  disturbance,  for  40  kHz  and  9i  =  80°.  x  is  the  direction  of  propagation 
and  y  is  the  direction  orthogonal  to  x,z  and  perpendicular  to  the  direction  of  propagation. 


gure  16) 


Figure  18.  Comparison  of  analytic  reflection  altitude  (Hr),  corrected  (Hrq),  and  altitude 

region  where  effective  altitude  should  be  selected  for  40  kHz,  0^  =  80°,  and  for 
the  range  away  from  a  high-altitude  burst  where  gamma  deposition  is  dominate. 


where 
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The  effective  reflection  region  is  still  lower  than  Hr^  for  much  of  the  disturbed 
region  shown  in  Figure  18.  However,  for  the  overall  set  of  computations  performed 
to  test  Equation  10,  it  performed  well. 


3.3  COMPUTATION  CRITERIA  FOR  ADJACENT 

VERTICAL  IONOSPHERIC  PROFILES 

It  has  been  long  known  that  the  reflection  coefficient  sensitivity  to 
changes  to  the  earth's  magnetic  field  characteristics  decreases  as: 


1.  The  ionosphere  becomes  more  depressed  due  to  increased  solar 
activity  (day  versus  night,  PCA  events)  or  to  nuclear 
disturbances . 

2.  The  distance  from  the  magnetic  equator  (increasing  magnetic 
dip  angle) . 

Correspondingly,  the  reflection  coefficient  sensitivity  to  changes  in  the  magnetic 
azimuth  also  decreases.  In  the  past,  WEDCOM  IV  used  a  broad  criteria  to  reduce  the 
number  of  computations  resulting  from  a  reduced  magnetic  field  anisotropy.  The 
criteria  relied  on  practical  experience  with  reflection  coefficient  computations 
where  it  was  known  that  the  magnetic  field  sensitivity  was  negligible.  The  dif¬ 
ficulty  with  this  criteria  was  the  broadness  of  it — that  it  permitted  repetition 
of  many  near  identical  reflection  coefficient  computations.  Therefore,  there  was  a 
need  to  develop  a  new  criterion  which  would  be  sensitive  to  magnetic  field  changes 
as  they  affect  the  reflection  coefficient  computation.  To  be  effective  this  new 
criterion  must  be  computationally  fast  as  was  the  old  criterion. 

A  function  used  in  the  reflection  coefficient  computation  is 


q 


,1/2 


♦i  + 


Xn 


(12) 


where  w  is  the  signal  frequency  (radians),  c  is  the  velocity  of  light,  and  <Jk  is  the 
ionospheric  incidence  angle  (assumed  here  to  be  80°) .  The  parameter  Xh  is  the  first 
diagonal  term  of  the  susceptibility  matrix  defined  as 


Xn  =  - 


X 

U 


£2y2\ 

-  Y2  / 


(13) 
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where 


i  -  yr 

X  =  (a)  /go)  ^ 

P 

U  =  1.0  j (v/w) 


Y  =  W  a)  /u)  (henrys) 
u  =  magnetic  gyro frequency  (s  *) 
ojp  =  plasma  frequency  (s  1 ) 

£  =  cos  0_  cos  d> 

D  ra 

0^  *  magnetic  dip  angle  (deg),  and 
<{>  =  magnetic  azimuth  (deg)  . 

3. 

For  an  isotropic  ionosphere  q  reduces  to  the  vertical  planewave  propagation  con¬ 
stant.  The  imaginary  part  of  q  is  thus  a  measure  of  the  signal  attenuation  at  the 
altitude  it  was  determined.  The  q's  computed  in  the  reflection  region  have  the 
sensitivity  to  the  magnetic  field  strength  important  to  the  reflection  coefficient 
computation . 

The  following  procedure  will  be  used  to  determine  whether  adjacent  loca¬ 
tions  along  the  propagation  path  should  have  separate  reflection  coefficient  com¬ 
putations  : 


1.  Determine  if  the  reflection  altitude,  H^,  has  changed.  If  so, 
compute  new  reflection  coefficients.  If  not, 

2.  Compute  q  for  both  vertical  profiles.  If  the  difference  in 
Im(q)  exceeds  6^  where 

JO. 03  for  Normal  Night 

6  = 

^  <0.005  for  all  other  cases 


(14) 


then  the  computations  for  both  locations  should  be  made. 

The  values  assigned  to  <S^  are  somewhat  arbitrary  and  must  be  considered 
preliminary.  They  were  obtained  by  comparing  the  value  of  q  with  the  reflection 
coefficient  behavior  for  the  cases  listed  in  Table  3.  Figures  19  through  21  are 
a  sample  of  the  data  generated.  Figure  19  illustrates  the  general  behavior  of  Im(q) 
for  magnetic  dip  and  azimuth  and  for  normal  and  a  moderately  disturbed  environment. 
Figures  20  and  21  are  sample  outputs  of  the  reflection  coefficient  sensitivity  to 
magnetic  field  effects  and  disturbed  environments.  Tables  of  representative  azimuth 
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Table  3.  Parameters  varied  in  q  and  ionospheric 
reflection  coefficient  calculations  to 
obtain  6  . 


PROFILE  Normal  Night  and  Normal  Day 

Spread  debris  night  =  10"^,  10 
10'10,  10"9 

DIP  ANGLE  30,  50,  70  degrees 

MAGNETIC  AZIMUTH  -85,  -65,  -45,  -25,  -5,  15,  35 
55,  75  degrees 

FREQUENCIES  27,  40,  100  kHz 


increments  (A)  were  created  where  A  represented  the  azimuth  change  where  the  magni¬ 
tude  of  „Rm  was  always  within  1  dB  change.  Note  the  sensitivity  of  A  with  dip  angle 
and  with  the  disturbed  environment.  The  proper  choice  of  6^  will  incorporate  this 
sensitivity.  The  <5^  values  in  Equation  14  are  not  stringent  and  will  not  have  the 
needed  sensitivity  for  some  cases  (especially  light  disturbed  cases).  A  more  con¬ 
servative  value  6  =  0.0001  will  assure  that  the  same  criterion  is  met  in  almost 

q 

all  cases.  6^  tends  to  require  a  larger  value  for  profiles  nearer  normal  night  and 
higher  frequencies  and  a  more  sophisticated  set  of  values  will  be  selected. 

3.4  INTERPOLATION  ALGORITHM 

In  the  past,  an  analytic  form  of  the  reflection  coefficient 

r i Cc i )  =  -  exp(aCi) 

C.  =  cos  0.  (15) 

ii 

has  been  used  (see,  for  example.  Reference  11  where  r^  represents  „RM,  iR1,  „RA 
or  ^R„,  and  0^  is  the  angle  of  incidence  at  the  ionosphere.  Only  a  single  reflection 
coefficient  computation  is  needed  to  determine  the  complex  valued  a.  This  form  was 
accurate  for  lower  VLF  and  for  sharp  gradient  isotropic  ionospheres.  The  isotropic 
formulation  in  WEDCOM  (Reference  1)  improved  this  form,  ie. 


ri(C.)  =  -  exp(aC.)  expjjo^C.J  (: 

where  the  real  valued  oij  better  accounted  for  less  sharply  reflecting  gradients. 
Neither  of  these  forms  did  well  in  representing  highly  anisotropic  profiles,  how¬ 
ever,  especially  where  the  reflection  gradients  were  not  simple. 
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Figure  20.  Magnitude  of  „R„  sensitivity  to  magnetic  field  for  normal  nighttime  and 


Figure  21.  Magnitude  of  |,R„  sensitivity  to  magnetic  field  for  moderate  spread  debris 
case  for  frequency  =  100  kHz. 


Two  recent  efforts  have  provided  more  information  on  representing  the 
reflection  coefficient  in  an  analytic  form  and/or  interpolating  between  values.  The 
work  described  in  Reference  2  utilized  an  analytic  form  of  the  reflection  coefficient, 
but  this  form  could  only  be  used  every  1/2  degree  (or  smaller)  in  6^  to  retain  the 
needed  accuracy.  Fortunately,  this  fine  increment  was  needed  only  for  the  MF  fre¬ 
quencies  being  analyzed,  and  it  was  shown  that  lower  VLF/LF  frequencies  did  not  re¬ 
quire  as  small  an  increment.  Figure  22  provides  examples  of  the  0^  sensitivity  for 
a  normal  nighttime  ionosphere  where  computations  were  performed  every  1  degree.  For 
the  MF  case  (350  kHz)  the  1-degree  increments  used  do  not  provide  the  needed  sensi¬ 
tivity  for  ARA.  For  the  LF  case  (100  kHz)  the  1-degree  increment  is  satisfactory. 

The  structure  noted  in  these  figures  become  less  apparent  as  the  frequency  is  re¬ 
duced  and  as  the  ionosphere  becomes  more  depressed.  It  is  interesting  to  note  the 
relative  values  of  the  diagonal  and  cross  terms  of  the  reflection  coefficient  matrix. 
It  is  difficult  to  generalize  the  characteristics  for  anisotropic  profiles  such  as 
this  one.  In  any  case  it  was  determined  that  an  analytic  form 

r.(C.)  =  exp(A+ BC.)  (17) 

would  well  represent  the  individual  anisotropic  reflection  coefficients  if  the  0^ 
increment  was  small  enough. 

The  second  work  which  required  an  analytic  form  of  the  reflection  coef¬ 
ficient  is  described  in  Reference  3.  A  very  large  data  base  of  reflection  coef¬ 
ficients  was  generated  parametrically  with  frequency  (27  through  100  kHz),  0^  (75 
through  82  degrees),  magnetic  dip  (30,  50,  70  degrees),  and  azimuth  (-75  through  85 
degrees),  day  and  night,  and  spread  debris  and  gamma  source  nuclear  environments. 

The  Equation  17  formulation  was  used  with  a  least-means-squared  measure  (as  a  function 
of  0^)  to  compute  representative  complex  valued  A  and  B.  With  few  exceptions,  the 
analytic  form  well  represented  the  reflection  coefficient,  especially  for  the  lower 
frequencies  and  for  the  more  disturbed  cases.  The  appropriateness  of  this  form  is 
seen  in  Figures  23  through  27.  Here  samples  of  the  general  linear  behavior  of  „R„ 
and  aRa  are  demonstrated  for  normal  and  disturbed  ionospheric  profiles,  frequencies 
of  27  and  40  kHz,  and  dip  angles  of  30  and  70  degrees.  The  phase  behavior  illus¬ 
trated  in  Figure  24  is  typical  (note  2n  discontinuity).  That  is,  it  has  a  very 
linear  behavior  with  cos  0^.  Note  the  change  in  the  sign  of  the  slope  between  |(R n 
and  aRa.  This  was  often  observed  for  the  anisotropic  profiles.  Figure  27  shows  how 
the  anisotropic  sensitivity  disappears  as  the  ionosphere  becomes  depressed. 


qure  22.  Reflection  coefficients  for  normal  nighttime  as  a  function  of  0. 


Figure  23.  Magnitude  of  WR„  and  XR  for  a  normal  night,  dip  =  70 
and  frequency  =  27  kHz. 
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Figure  24.  Phase  of  HR„  and  ^  for  normal  night,  dip  =  70 
and  frequency  =  27  kHz. 
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and  frequency  =  27  kHz . 
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Figure  26.  Magnitude  of  „R„  and  XRX  for  normal  night,  dip 
and  frequency  =  40  kHz. 


The  analytic  form  of  the  reflection  coefficient,  liquation  17,  will  he  used 
in  the  revised  WEDCOM.  This  permits  a  more  efficient  EE  and  VI. E  model  simulation 
with  little  sacrifice  in  accuracy.  Section  3.3  describes  a  procedure  to  eliminate 
unnecessary  reflection  coefficient  computations  if  the  ionospheric  profile  and 
anisotropy  has  changed  negligibly.  Reference  3  describes  a  reflection  coefficient 
scaling  procedure  which  is  slightly  modified  here  to  account  for  anisotropic  changes. 
When  this  scaling  is  used  in  conjunction  with  the  analytic  formulation  just  described, 
a  first-order  correction  is  made  to  the  reflection  coefficient  for  a  slightly  mod¬ 
ified  ionospheric  profile.  The  simplicity  of  this  scaling  makes  it  desirable.  How¬ 
ever,  as  pointed  out  in  Reference  3,  this  scaling  works  well  only  when  the  ionos¬ 
pheric  profile  to  be  scaled  is  very  similar  to  the  one  for  which  the  reflection  coef¬ 
ficient  was  computed.  This  requirement  should  be  satisfied  by  the  Section  3.3  pro¬ 
cedure.  The  susceptibility  term.  Equation  13,  is  assumed  to  behave  exponentially  in 
altitude  with  a  scale-height  H  in  the  reflection  region.  The  reflection  coefficient 
for  the  new  profile  is  scaled  from  the  complex  A  and  B  of  the  present  profile  from 


ri(Ci>  =  «p(a  ♦  B  C.) 


(  18) 


where 


H  ^  =  susceptibility  scale-height  for  new  profile 

H^p  =  susceptibility  scale-height  for  present  profile. 

This  scaling  procedure  is  emperical,  but  has  its  foundation  in  an  analytic  treatment 
of  exponentially  varying  isotropic  ionospheres  described  in  Reference  12. 

The  analytic  formulation  (Equation  13)  will  not  significantly  reduce  the 
computation  time  if  a  large  number  of  iterative  integration  reflection  coefficient 
computations  (Section  3.1)  are  needed  to  determine  the  complex-valued  A  and  B  con¬ 
stants.  Furthermore,  the  required  number  of  computations  vary  as  a  function  of  the 
ionospheric  environment  and  frequency.  For  this  analysis  the  A  and  B  constants  were 
determined  by  first  computing  the  reflection  coefficients  at  two  incidence  angles 
and  then  solving  for  A  and  B.  This  procedure  was  noted  to  have  several  limitations, 
especially  in  accounting  accurately  for  the  2tt  ambiguities  as  noted  in  Figure  24. 
Also,  the  region  of  accuracy  is  unknown.  These  limitations  were  not  important  in 
the  EE  analysis,  but  became  critical  in  the  VEE  study. 

A  more  accurate  procedure  is  to  first  determine  the  reflection  coefficient 
(R)  and  its  derivative  (dR/dC.)  at  a  specified  incidence  angle  ((!.).  Then  by  noting 


dR/dC. 


(19) 


after  which  A  is  easily  evaluated,  the  phase  ambiguity  is  eliminated.  By  computing 
B  in  this  manner  for  one  or  two  more  Ch  ,  the  A's  and  B's  can  be  compared  to  note  the 
sensitivity.  If  needed,  intermediate  Ch  computations  are  then  made. 
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SECTION  4 

VLF  MODEL  IMPROVEMENTS 


4.1  INTRODUCTION 

The  WEDCOM  VLF  model  underwent  a  major  revision  from  WEDCOM  III  to  the 
present  WEDCOM  IV  version.  The  most  important  improvement  in  WEDCOM  IV  was  the 
addition  of  an  anisotropic  ionosphere  capability,  poorly  approximated  in  the  earlier 
model.  This  required  a  much  larger  computation  time,  however.  For  example,  the 
older  version  computed  only  two  (or  three)  ionospheric  reflection  coefficients  for  a 
normal  nighttime  case  whereas  the  present  version  typically  computes  36  to  50  values. 
The  older  version  required  only  two  or  three  values  due  to: 

1.  Use  of  an  ionospheric  reflection  coefficient  interpolation 
algorithm, 

2.  Fewer  reflection  coefficient  segments. 

The  penalty  of  the  older  version  was  that  the  VLF  predictions  were  not  accurate  for 
nighttime  conditions.  However,  the  older  version  provided  accurate  predictions  for 
highly  nuclear  disturbed  environments,  and  was  not  too  much  in  error  for  normal 
daytime  conditions. 

In  addition  to  the  reflection  coefficient  computation,  the  WEDCOM  IV  pro¬ 
cedure  that  determines  the  eigenangles  is  more  sophisticated  and  also  required  more 
computation  time  than  the  earlier  procedure.  As  the  WEDCOM  IV  VLF  model  has  been 
exercised,  it  has  become  apparent  that  the  increased  computation  time  is  excessive 
for  the  objectives  of  WEDCOM,  and  at  least  a  portion  of  this  increased  time  is  for 
a  computational  accuracy  not  needed  for  WEDCOM  predictions.  The  alternative  model 
modifications  presented  below  suggest  reductions  in  the  computation  time  while  re¬ 
taining  the' needed  accuracy  and  model  capability. 

Since  WEDCOM  IV,  the  VLF  propagation  model  was  revised  for  use  in  studies 
for  the  Naval  Electronics  System  Command  (Reference  2).  The  revision  includes  the 
addition  of  a  more  accurate  computation  for  highly  elevated  antennas  and  a  redefini¬ 
tion  of  the  eigenangle  search  region.  The  first  improvement  is  not  necessary  for  the 
purposes  of  WEDCOM  and  will  not  be  included.  The  second  improvement  which  allows  a 
reduction  in  the  number  of  reflection  coefficients  needed  is  discussed  in  Section  4.2. 
Additional  improvements  were  also  made  that  provide  the  model  increased  capability. 


For  example,  an  improved  mode  convergence  computation  has  been  developed,  and  a 
variable  path  increment  logic  was  added  to  reduce  the  number  of  computations. 

There  are  two  areas  where  large  computation  time  requirements  are  needed 
for  the  VLF  model: 

1.  The  ionospheric  reflection  coefficient  computation,  and 

2.  The  search  for  the  eigenangles. 

The  reflection  coefficient  computation  was  addressed  in  Section  3,  where  an  analytic 
formulation  to  reduce  the  time  requirement  was  described.  An  improved  criteria  for 
when  the  computations  are  to  be  performed  (also  described  in  Section  3)  will  also  re¬ 
duce  the  computation  time.  The  second  area,  the  eigenangle  search,  becomes  increas¬ 
ingly  important  as  the  number  of  eigenangles  increase.  Alternative  procedures  for 
eigenangle  determinations  are  described  in  Sections  4.3  and  4.4. 

The  new  path  increment  procedure  outlined  in  Section  2.2  for  the  LF  model 
will  be  utilized  in  the  VLF  model  to  determine  the  eigenangles. 

4.2  PROCEDURE  TO  REDUCE  THE  NUMBER  OF 

REFLECTION  COEFFICIENT  COMPUTATIONS 

WEDCOM  IV  used  two  techniques  for  determining  the  eigenvalues,  8^  (see 
Reference  1) .  The  search  technique  over  the  0-plane  is  the  major  user  of  computer 
time  and  is  the  one  addressed  here.  It  is  an  adaptation  of  the  MODE  SRCH  algorithm 
(Reference  13).  Figure  28  illustrates  the  WEDCOM  IV  method  of  search  "boxes"  where 
an  ionospheric  reflection  coefficient  computation  is  made  at  each  corner  (ie,  for  N 
boxes,  there  are  4N  computations).  This  overlapping  of  the  boxes  assured  that  eigen¬ 
values  falling  on  the  box  edges  would  not  be  missed.  Figure  28  illustrates  the 
modified  procedure  which  reduces  the  number  of  computations  at  least  by  25  percent, 
but  sometimes  approaching  50  percent.  In  utilizing  this  procedure  there  were  a  few 
occassions  when  the  eigenvalue  did  fall  on  a  box  edge  and  was  not  computed.  For 
these  few  cases  a  special  computation  was  made  to  determine  the  eigenangle,  which 
worked  satisfactorily  for  those  exercises.  For  WEDCOM,  a  user  interaction  such  as 
this  is  not  practical  and  a  simple  procedure  has  been  outlined  (but  not  yet  imple¬ 
mented)  to  eliminate  this  difficulty.  In  essence  the  search  boxes  will  be  slightly 
enlarged  while  retaining  the  same  reflection  coefficient  locations.  The  duplication 
of  an  eigenvalue  by  this  procedure  will  be  automatically  identified  by  the  present 
WEDCOM  IV  algorithm. 

An  additional  computation  time-saving  technique  was  implemented  into  the 
VLF  model  for  the  NESC  studies.  For  cases  where  ionospheric  conditions  are  identical 
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Figure  28.  Comparison  of  new  and  WEDCOM  IV  procedures  for 

determining  0-plane  reflection  coefficient  computations. 


for  adjacent  path  segments,  but  where  the  ground  conductivity  does  change,  the 
ionospheric  reflection  coefficients  are  not  recomputed  as  was  done  in  WEDCOM  IV. 

The  degree  of  computational  time  savings  for  this  procedure  is  obviously  case 
dependent . 

4.3  PRECOMPUTED  EIGENANGLES 

A  recent  effort  (Reference  3)  created  an  eigenangle  data  base  for  normal 
and  disturbed  environments  from  which  needed  eigenangles  were  obtained  through 
interpolation  and  scaling  algorithms.  From  this  procedure,  extremely  fast  evalua¬ 
tions  (even  compared  to  the  WEDCOM  III  model)  are  possible.  The  amount  of  subroutine 
support  is  a  very  small  fraction  of  WEDCOM  IV.  The  price  for  this  efficiency,  how¬ 
ever,  is  an  immense  data  file  and  a  computational  accuracy  below  the  WEDCOM  standards 

The  most  difficult  and  time  consuming  eigenangle  evaluations  occur  near 
normal  nighttime.  A  proposed  technique  is  to  have  a  precomputed  data  file  only  for 
this  case.  Alternative  methods  would  be  used  for  all  other  cases.  The  data  file 


would  still  be  very  large  but  would  be  manageable.  The  data  base  could  be  generated 
by  two  procedures: 

1.  Precomputed  and  furnished  with  the  WEDCOM  Fortran  tape. 

2.  A  separate  computation  performed  by  the  user. 

Initially,  the  first  procedure  may  seem  the  simpler.  However,  experience  with  this 
procedure  has  indicated  frequent  difficulty  with  the  data  tape  and  an  inconvenience 
in  generating  a  new  tape  when  model  modifications  are  required.  The  second  pro¬ 
cedure  is  favored  at  the  present,  where  a  separate  driver  would  be  furnished  to  the 
user.  This  would  also  allow  the  data  to  better  conform  to  the  user  needs. 


4.4  EIGENANGLE  DETERMINATION  USING  ANALYTIC 

REFLECTION  COEFFICIENT  FORMULATION  AND 
INITIAL  EIGENANGLE  ESTIMATES 

The  mode  search  model  used  in  WEDCOM  IV  to  find  an  initial  value  for  the 
eigenangle  solution  requires  significant  calculation  time,  particularly  for  lightly- 
disturbed,  nighttime  conditions.  In  WEDCOM  III  analytic  r  .proximations  were  used 
for  the  initial  eigenangles.  A  preliminary  study  has  been  made  to  see  if  a  similar 
procedure  (using  the  analytic  form  for  the  reflection  coefficient  described  in  Sec¬ 
tion  3)  can  be  used  to  replace  the  mode  search  model.  A  summary  of  the  eigenangle 
approximations  follow.  Each  approximation  applies  to  one  of  three  easily  identifi¬ 
able  eigenvalue  types.  The  eigenangle,  0N>  is  defined  at  altitude  HR. 

4.4.1  Quasi -FI at  Earth  Modes 


=  COS  0, 


where 


2  .OkH 


lR  -  j  B  [l .  0 


+  4.0  HR(kHR)V(aRj 


'jA.0(kHR)2  +  B2! 


(2  .OttN  -  j  A 


(vertical  polarization) 


'tt(2  .ON  +  1 .0)  -  jA  (horizontal  polarization) 

N  is  the  mode  number,  and  k,  a,  HR,  A,  and  B  were  defined  in  Section  3. 

4.4.2  Earth-Grazing  Modes 

c;  =  [cf  +  xl1/2  =  cos(0J 


where 


■e?) 


Cl  -  \h  -  (2.0ka/3.0)X^  -  jBXa|  /  jk^  ♦  j 


tt(12 .ON  +  1 .0)/6 .0  -  jA  (vertical  polarization) 


tt(12.0N  + 6.0)/6.0  -  jA  (horizontal  polarization) 


4.4.3  Whispering-Gallery  Modes 


where 


’  '  =  Y  )  1 


:i-[c;*xa]‘ 


0  -  xJj  -  k. 


cos(0N) 


k  X  +  j 
a  a  J  X 


-  JBxJ/jk  x  .  j 


X2  =  Rn  -  (2.0ka/3.0)X-  -  jBXaj/jkaXa 


Rn  =  tt(4  .ON  +  1 ) / 2 . 0  -  jA  .  (29) 

Note  that  the  values  of  A  and  B  will  differ  dependent  on  the  polarization.  Infinite 
earth  conductivity  is  assumed  in  the  above  approximations.  The  best  eigenvalue  ap¬ 
proximation  is  selected  based  on  its  type.  Then  an  iteration  procedure,  using  the 
isotropic  definition  of  the  model  solution,  is  used  to  improve  the  approximation 
sensitivity  to  the  actual  earth  conductivity.  The  resulting  eigenvalue  approxima¬ 
tion  is  iterated  to  a  final  exact  solution  using  the  anisotropic  model  solution. 

In  WEDCOM  IV  the  reflection  coefficients  for  an  isotropic  ionosphere  are 
used  (earth-curvature  and  magnetic  field  effects  are  ignored).  In  the  new  model  the 
approximate  eigenangle  formulations  will  retain  the  isotropic  ionosphere  assumption, 
but  will  use  anisotropic  reflection  coefficients.  Figure  29  is  reproduced  from 
Reference  1  and  demonstrates  how  the  approximate  isotropic  procedure  is  usually  ac¬ 
curate  for  daytime  and  strong  nuclear  depressed  ionospheres.  As  a  normal  nighttime 
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Figure  29.  Loci  of  eigenangles  for  the  approximate  isotropic  and 
exact  anisotropic  procedures. 

is  approached  the  first  TE  and  TM  modes  converge.  Figure  29  also  shows  why  the 
older  analytic  formulation  had  a  built-in  distortion.  The  ionosphere  reflection 
coefficients  used  as  data  points  for  the  isotropic  analytic  formulation  are  computed 
along  the  real  0  axis  (this  is  also  the  case  for  the  LF  model).  Especially  for  the 
lower-ordered  modes  the  eigenangles  are  not  near  the  real-axis  (and  thus  not  near  the 
data  points) .  For  the  model  effort  described  here,  complex  valued  0  were  used  as 
data  points.  The  imaginary  part  of  the  0  value  was  selected  to  correspond  to 
1  dB/Mm  absorption. 

The  new  formulation  has  had  limited  testing.  A  normal  daytime  case  was 
executed  and  resulted  in  significant  computation  efficiency  compared  to  WEDCOM  IV, 
and  had  surprisingly  good  accuracy.  A  normal  nighttime  case  (the  more  interesting 
and  difficult  case)  was  then  exercised.  For  this  case  the  approximate  eigenangle 
did  not  always  provide  an  accurate  estimate  for  proper  convergence  in  the  anisotropic 
solution.  Several  factors  contributed  to  this: 

1.  One  or  both  first-order  modes  sometimes  switch  polarization. 

2.  The  analytic  formulation  of  the  ionosphere  were  not  correct. 


3.  The  isotropic  mode  solution  did  not  sufficiently  approximate  the 
anisotropic  solution. 

At  this  time  the  second  item  appears  to  be  the  most  dominant  problem.  The  reflection 
coefficient  phase  values  were  not  always  correct  leading  to  erroneous  approximations. 
An  improved  procedure  to  correct  this  has  been  outlined.  An  alternate  procedure  has 
also  been  outlined  to  address  the  other  two  difficulties.  Briefly,  this  procedure 
"iterates"  from  the  approximate  value  (with  its  model  equation)  to  the  actual  eigen¬ 
value  (with  the  polarizations  and  anisotropy  of  the  actual  model  equation).  This 
requires  additional  calculation  time,  however,  and  must  be  evaluated  in  comparison 
to  the  precomputed  eigenvalue  procedure  (Section  4.3). 


SECTION  5 

ELF  PROPAGATION  MODEL 


5.1  GENERAL 

The  ELF  propagation  model  currently  in  the  WEDCOM  code  is  essentially  the 
same  as  the  Naval  Ocean  Systems  Center  (NOSC)  model  described  by  Pappert  and  Moler 
(Reference  14).  In  this  model  mode  solutions  for  propagation  in  the  earth- ionosphere 
waveguide  are  found.  Waveguide  variability  along  the  great  circle  propagation  path 
is  accounted  for  (WKB  approximation)  but  variability  normal  to  the  great  circle  path 
is  neglected.  Allowance  is  made  for  propagation  over  both  the  short  and  long  great 
circle  paths. 

The  iterative  Newton-Raphson  method  is  used  to  find  the  mode  solution. 

The  solution  is  started  with  a  guess  for  the  eigenangle.  Only  one  mode  propagates 
at  ELF  and  it  has  been  found  that  almost  any  reasonable  initial  guess  for  the  eigen¬ 
angle  will  result  in  convergence.  A  significant  part  of  the  mode  solution  in  terms 
of  computation  time  is  the  calculation  of  the  reflection  coefficient  matrix.  The 
matrix  is  found  from  a  numerical  integration  of  differential  equations  defined  by 
Budden  (Reference  15) .  For  ELF  the  integration  must  start  at  fairly  high  altitudes 
(about  200  km  for  ambient  nighttime  conditions)  and  this  increases  the  computation 
time  over  that  for  VLF  and  LF  calculations  where  the  integration  can  be  started 
below  100  km. 

Recently,  Booker  (Reference  16)  has  proposed  an  approximate  model  for  ELF 
propagation  that  combines  theoretical  work  of  Booker  and  Lefeuvre  (Reference  17) 
and  Griefinger  and  Griefinger  (References  18  and  19).  The  model  is  based  on  dividing 
the  atmosphere  into  regions  where  the  local  vertical  gradient  is  large  in  comparison 
to  the  local  wavelength  and  regions  where  it  is  small  in  comparison  to  the  local  wave¬ 
length.  In  regions  where  the  gradient  is  small  the  phase  integral  method  is  used  to 
determine  propagation.  Regions  where  the  gradient  is  high  are  replaced  with  a  dis¬ 
continuity  equal  to  the  difference  in  refractive  indices  at  the  top  and  bottom  of 
the  region.  Up  to  five  reflection  regions  can  be  defined  to  account  for  reflections 
from  the  D  and  E  regions.  Reflections  from  the  discontinuities  are  combined  using 
the  phase  change  between  the  reflection  regions  (phase  change  in  regions  where  phase 
integral  method  is  applicable).  As  described  by  Booker  phase  and  group  propagation 
below  the  bottom  of  the  ionosphere  is  practically  horizontal  while  phase  propagation 


above  the  bottom  of  the  ionosphere  is  practically  vertical.  The  bottom  of  the  ionos¬ 
phere  is  defined  by  a  complex  height.  The  positive  imaginary  part  is  a  measure  of 
the  absorption  experienced  in  the  lowest  layer  of  the  ionosphere. 

The  approximate  model  requires  between  1/10  and  1/5  the  computation  time 
required  by  the  detailed  numerical  model.  This  significant  reduction  in  computation 
time  will  allow  useful  system  analyses  for  conditions  that  are  too  expensive  with  the 
current  WEDCOM  code . 

5.2  BOOKER  MODEL 

A  detailed  description  of  the  model  developed  by  Booker  is  given  in  Refer¬ 
ence  16.  The  following  summarizes  the  computational  steps  and  principal  equations. 
Calculations  are  made  for  both  the  ordinary  (0)  and  extraordinary  (X)  waves.  Sub¬ 
scripts  identify  quantities  for  the  0  and  X  waves  when  results  for  both  waves  are 
combined . 


5.2.1  Evaluate  the  Heights  of  Reflection 

The  reflection  heights  are  determined  as  the  heights  where 


2tt  ]  n  |  h 


where 


X  =  free-space  wavelength 

n  =  complex  index  of  refraction  for  vertical  propagation 
h  =  1 1/n2  -  1  d/dz  (n2  -  1) I" 1 

and  y  is  a  quantity  that  is  a  function  of  the  refractive  index  phase  angle  (given 
graphically  in  Reference  16).  For  most  conditions  y  can  be  taken  as  2  for  the  0 
wave  and  2.57  for  the  X  wave.  Up  to  five  reflection  heights  can  be  computed  for  the 
0  wave  (h1Q>  h2Q,  h3Q  h4Q,  and  h50)  and  the  X  wave  (hlx>  h2x,  h3X,  h4X>  and  h5X) . 

The  complex  index  of  refraction  indices  are  found  from 
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a>Nm+ 


=  5.98  x  10  N 


<4n- 


=  5.98  x  10  N 


“Me 

<*Vlm+ 

“Mm- 


=  -1.76  x  10  B 


=  3.3  x  10  B 


=  -3.3  x  10  B 


-3, 


=  electron  density  (cm  ) 


-3. 


=  positive  ion  density  (cm  ) 


-3. 


=  negative  ion  density  (cm  ) 


B  =  magnetic  field  strength  (gauss) 


-1. 


=  electron-neutral  collision  frequency  (s  ) 


v. 
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-1. 


=  ion-neutral  collision  frequency  (s  ) . 


Positive  and  negative  atomic  ions  can  also  be  included  in  the  susceptibility  terms 
but  are  neglected  in  the  ELF  calculations.  For  computational  purposes,  Equation  31 
is  rewritten  as  (Reference  20) 

2^ 
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(35 


where 


2  2 

Y  =  -2(1 +<L)(1 +<T)  +  (1  +  <T)  (kl  -  <T)  -  <H)  cos  I 


S  =  y  -  4ae 


2 

a  =  1  +  -  (<L  -  k,j.)  cos  I 

e  =  (1  +  K,)  ((1  + <T) 2  + 

5.2.2  Calculate  the  Complex  Index  of  Refraction 
Indices  at  the  Reflection  Heights 

The  equations  given  in  Section  5.2.2  are  used  to  compute  the  indices.  Th 


The  complex  phase  changes  are  given  by: 


<&12  =  (2tt/A)  ndz  , 


$23  =  C27T/A)  (z3  -  z2) 


$34  =  (2tt/A)  /  H  ndz 


5.2.5 


*45  =  (2tt/A)(z5  -  z4)  . 


Calculate  the  Equivalent  Complex 
Height  of  the  Bottom  of  the  Ionosphere 

The  equivalent  height  of  the  bottom  of  the  ionosphere  is  computed  from 
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J0 


5.2.6 


Calculate  the  Phase  Velocity 
and  Attenuation  Rate 


The  phase  velocity  and  attenuation  rate  are  computed  from 
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where  R(S)  and  I(S)  are  the  real  and  imaginary  parts  of  S  and  S  is  given  by 
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5.2.7  Interpolation  and  Cal¬ 
culation  of  Derivatives 

In  computing  the  reflection  heights  it  is  necessary  to  determine  the  mag¬ 
nitude  and  height-derivatives  of  electron  and  ion  densities,  collision  frequencies, 
and  magnetic  field  strength.  As  pointed  out  by  Booker,  confusion  results  if  there 
are  discontinuities  in  the  magnitude  or  derivatives  of  any  of  the  quantities.  In 
choosing  the  electron  and  ion  density  profiles  Booker  used  a  method  he  developed 
(Reference  21)  of  fitting  the  profile  with  exponential  fits  between  selected  alti¬ 
tudes  (transition  altitudes)  and  then  using  smoothing  functions  at  the  transition 
altitudes  that  smooth  the  profile  over  a  prescribed  scale  length.  The  relation  for 
electron  and  ion  density  presented  by  Booker  is 


in 

£nN(z)  =  £nN  +  A  (z  -  z  )  +  N ^  (A  )  -  A  ) 
v  ’  r  0,1  /  j  v  n,n+l  n-l,n 


f(z  -  z  ,B  )  -  f(z  -  z  ,B  ' 
I  ^  r’  nJ  1  r  n’  nt 


where 


z  =  altitude 

zn  =  transition  altitudes  used  to  describe  electron  density 

A  ,  =  slope  of  £nN(z)  versus  z  for  z  <  z  <  z  , 
n,n+l  r  v  j  n  n+1 

B^  =  reciprocal  of  smoothing  scale  used  at  transition  altitudes 


f(z,B)  = 


B  1  + exp(Bz) j  zB  < 


zB  >  100 


N  =  N(z  ) 

r  r 

z  =  reference  altitude, 

r 

For  use  in  the  WliDCOM  code  where  electron  and  ion  densities  profiles  are 
computed  as  a  function  of  nuclear  radiation  it  is  convenient  to  use  Booker's  formula 
tion  as  an  interpolation  procedure.  This  can  be  done  by  choosing  the  smoothing 
scales  so  that  the  effects  of  fitting  outside  the  interpolation  region  are  negligibl 
within  the  interpolation  region.  A  4-point  interpolation  is  used  by  choosing  the 
smoothing  scale  so  that  the  profile  within  each  altitude  interval  is  determined  by 
the  exponential  fits  for  the  interval  and  for  the  adjacent  altitude  intervals.  Thus 
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N(z)  between  z 2  and  z3  is  given  by 
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The  derivative  of  N(z)  can  be  obtained  from 
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5.3  EVALUATION  OF  BOOKER  MODEL 

5.3.1  Comparison  with  Results  Presented 

by  Booker  and  Behroozi-Toasi 

Selected  comparisons  were  made  with  model  results  presented  in  Reference 
16  to  verify  that  the  implementation  of  the  Booker  model  was  correct.  Figures  30 
through  33  show  comparisons  of  results  for  phase  velocity  and  attenuation  rate  for 
quiet  ionosphere  profiles.  The  small  differences  are  due  to  the  different  methods 
used  for  interpolation  and  to  compute  derivatives.  Comparison  of  the  results  at  75 
degrees  north  latitude  with  detailed  calculations  by  NOSC  are  given  in  Reference  16. 
The  results  are  in  excellent  agreement. 

5.3.2  Comparison  with  Selected  NOSC  Results 

Calculations  with  the  Booker  model  were  also  made  for  selected  profiles 
where  detailed  NOSC  results  were  available.  The  profiles  selected  include  day  and 
night  quiet  ionosphere  profiles  currently  used  in  the  WEDCOM  IV  code,  day  and  night 
quiet  ionosphere  profiles  used  by  NOSC  in  Reference  14,  and  a  profile  representative 
of  moderate  nuclear  ionization.  Tables  4  through  8  show  the  electron  and  positive 
ion  densities  and  collision  frequencies  used. 

Tables  9  and  10  show  a  comparison  of  the  phase  velocity  and  attenuation 
rate  obtained  with  the  Booker  model  and  with  the  detailed  NOSC  model.  Calculations 
were  made  at  70-degrees  magnetic  north  latitude  (80-degrees  magnetic  dip  angle). 
There  is  good  agrement  between  the  models  except  for  nighttime  conditions  where  the 
Booker  model  overestimates  the  attenuation  rate.  The  cause  of  this  difference  is 
still  under  investigation. 

5.4  CALCULATION  OF  EXCITATION 
AND  HEIGHT-GAIN  FACTORS 

In  order  to  compute  electric  and  magnetic  field  strengths  it  is  necessary 
to  specify  excitation  and  height-gain  factors.  The  height-gain  factors  for  a  ver¬ 
tical  dipole  launch,  end-on  launch,  and  broadside  launch  from  a  horizontal  dipole 
are  (Reference  14) 
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RATIO  OF  VELOCITY  OF  LIGHT  TO  VELOCITY  OF  PHASE  PROPAGATION 


LATITUDE 

Chapter  4 

Figure  30.  Illustrating  the  variation  of  the  velocity  of  phase 
propagation  with  latitude  and  frequency  under  daytime 
conditions.  Solid  curves  reproduced  from  Reference  16. 
Symbols  indicate  values  obtained  with  Booker's  model 
at  Tempo. 
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RATIO  OR  VELGCfTY  OR  UQMT  TO  VElOCfTV 


Figure  32.  Illustrating  the  variation  of  the  velocity  of  phase 
propagation  with  latitude  and  frequency  under  night 
time  conditions.  Solid  curves  reproduced  from  Ref¬ 
erence  16.  Symbols  indicate  values  obtained  with 
Booker's  model  at  Tempo. 


Table  4.  Ionization  and  collision  frequencies  for 
WEDCOM  quiet  nighttime  ionosphere. 
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Table  5.  Ionization  and  collision  frequencies  for  NOSC  quiet 
nighttime  ionosphere. 
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2 , 60E ♦OP 
1 ,29Et09 
6,37Ef08 
3. 15E+08 
1 .55E+08 
7,70Ef07 
3.80E+07 
1 ,90Et07 
9.30E+06 
4.60E+06 
2.30E^06 
1 ,l0Et06 
5.50Eto5 
2, 70Et05 
1 ,34Et05 
6.64E+04 
3.28E+04 
1 .60E  +  04 
8.00E+03 
4.00E+03 
2.00E+03 
1  ,OOE^03 
3.00E+02 
1 ,50E^02 
8,00Et0l 
1 . 80E ♦ 0  1 
6 , 00E ♦OO 
2,00£t00 
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Table  6.  Ionization  and  collision  frequencies  for  WEDCOM 
quiet  daytime  ionosphere. 


z 

Ne  (cm'3) 

N+  (cm-3) 

Ve  (s'1) 

(s_1) 

ALT(KM) 

ENE(CM-3) 

EnP ( CH»3  j 

XNUE  C  S- 1 ) 

XNUI (S*i ) 

0,00 

1.71E-07 

3.51E+02 

1.53EM1 

6.00E+09 

5.00 

2.00E-06 

1.26Et03 

8.22EM0 

3.39£to9 

10,00 

1.37E-05 

2.44E+03 

4 , 1 1 E  + 1 0 

t  ,80Ef09 

15,00 

8,26E»05 

3.43E+03 

1 . 9 1  E  ♦  t  0 

8.62E+08 

20,00 

4.56E-04 

4, 05Ef03 

8.72E+09 

3.93E+08 

25,00 

2. 18E-03 

4 , 1 3E*03 

4.03E+09 

1 ,80Et08 

30.00 

8.80E-03 

3.69E+03 

1 .90E+09 

8.36E+07 

35,00 

3, 12E-03 

3.13E+03 

9, 19E  +  08 

3.95E+07 

40,00 

8. 14E-02 

2,45Et03 

4.62E+08 

1 ,94E*07 

45,00 

1 .94E-01 

1 ,82E*03 

2.41E+08 

9,87E*06 

50,00 

2.53Et00 

1 ,27E*03 

1 ,29Et08 

S.23E+06 

55,00 

1.79E+01 

6 . 1  1  E  ♦  0  2 

6.96E407 

3 , 26E+06 

60,00 

5.59E+01 

2.04E+02 

3.65E+07 

2,0lEt06 

65,00 

4.49E+02 

6.03E402 

1 ,86E*07 

1 .21E+06 

70.00 

1.76E+03 

1 .85E+03 

9.02E+06 

7 , 0 1 E*05 

75.00 

4,08e^03 

4.11E+03 

4.12E+06 

3.84Ef05 

80,00 

4,06E*O3 

4.86E+03 

1 .75E+06 

1 ,96E*05 

85,00 

5.64e^03 

5.64Ef03 

6.99E+05 

9.2lEt04 

90,00 

9.27E+03 

9.27E+03 

2.73E*05 

4.13E+04 

95,00 

4 , 55e  +  04 

4 , 55E ♦  04 

1 . 10E^05 

1 .87E+C4 

100.00 

8.39E+04 

e.39E+04 

4 .65E+04 

8.95E*03 

105,00 

1 ,2lEt05 

1 ,21E*05 

2.14E+04 

3.95£f03 

110.00 

1 • 54E ♦ 05 

1 .54E*05 

1 . 10E+04 

1 ,92E*03 

115.00 

1 ,85e  +  05 

1 ,85Et05 

6,22Ef03 

1 .02E+03 

120.00 

2 . 0 1  E*05 

2.01E+05 

3,87Et03 

5 . 89E+02 

130.00 

2.10E+05 

2. 10E+05 

1 .86E+03 

2.48Et02 

140.00 

2.30E+05 

2, 30E+05 

1 ,06Ef03 

1 ,33Et02 

160.00 

2.97E405 

2.97Et05 

4 » 32E  +  02 

5.27E+01 

180,00 

3.39Et05 

3.39E+05 

2, 04E+02 

2.57E+01 

200.00 

3.61E*05 

3,6lEt05 

l,04Et02 

1 . 4  0E ♦ 0  1 

Table  7.  Ionization  and  collision  frequencies  for  NOSC  quiet 
daytime  ionosphere. 


altikm) 

0,00 

5.00 

10,00 

15.00 

20.00 

25,00 

30,00 

35.00 

«0.OC 

05.00 

50.00 

55.00 

60,00 

65.00 

70.00 

75.00 

80,00 

85.00 

90,00 

95.00 

100,00 

105,00 

110.00 

115,00 

120,00 

130,00 

190,00 

160.00 

180,00 

200.00 


Ng  (cm-3) 

ENE(CM-3) 

1 .99E-06 

3.35E-06 

6.83E-06 

1 .55E-05 

3.54E-05 

8.45E-05 

1 .95E-04 

4.00E-04 

1 ,00E*03 

1 .00E-02 

1 .00E-01 

1 ,00Et00 

1 ,00E*01 

2.25Ef02 

7.00Et02 

1 ,50£f03 

3.00E+03 

6.00E+03 

1 ,20Et04 

3,00E*04 

6,O0EtO4 

1 , 00Et05 

1.20Et05 

1 .30E  +  05 

1 ,40Ef05 

1 ,60Et05 

2.00E+05 

4.00E+05 

6,O0E*05 

8,00Et05 


N+  (cm-3) 

ENP (CM*3) 
5 , 50E*03 
4.00E+03 
3.50E+03 
3.30E+03 
3,40E*03 
3.70Et03 
3.30E+03 
3.00E+03 
2.50E+03 
2.00E*03 
1 , 50Et03 
1 ,00E*03 
8.00E+02 
5.00Ef02 
9 . 00E ♦ 02 
1 ,50Ei03 
3.00E*03 
6.00E+03 
1 ,20E*04 
3.00E+04 
6.00E+04 
1.00E+05 
1 .20E  +  05 
1 ,30Et05 
1  ,«0Et05 
1 .60E  +  05 
2.00Et05 
4.00E+05 
6.00E+05 
8,00E*05 


ve  (s'1) 

XNUE(S«1  ) 
4 , 30Et  1 1 
1 ,90Etl  1 
8.50EM0 
3,77EtlO 
1 .68E  +  1  0 
7 , 45E ♦ 09 
3.30E+09 
1.47E+09 
6.50E+08 
2.70E+08 
1 .30E+08 
5.70E+07 
2.60Et07 
1 . l3Et07 
5.00E+06 
2.20E+06 
9.96E+05 
4.00E+05 
1 .97E+05 
8.70E+04 
3.90E+04 
2.60E+04 
1 .80E+04 
1 .04E+04 
1 ,00E*04 
4.00E+03 
2,50E>03 
6.00E+02 
1 ,70E*02 
4.50E+01 


Vi  (s'1) 

XNUI(S-l) 
2.14E+10 
1 . 06E ♦  l  0 
5.20E+09 
2.58E+09 
1 ,27Et09 
6.30E*08 
3. 10E*08 
1 ,54Et08 
7.60E+07 
3.80E»07 

1 , B6E*07 
9 , 20E+06 
4.fe0EtO6 
2 ,20E*06 
1 . 10E+06 
5.40E+05 

2, b8Et05 
1 ,33Et05 
6,56Et04 
3.20E+04 
1 .60E+04 
8.00E+03 
4.00E+03 
2.00E+03 
6S00E*02 
3, 00E+02 
1 .60E  +  02 
3 . 60E*0 1 
1 .20E  +  01 
4 .  OOE^OO 


Table  8.  Ionization  and  collision  frequencies  for 
moderate  nuclear  ionization. 


Ne  (cm-3) 


N+  (cm"  ) 


ve  (s_i) 


V.  (s'1) 


ALT (KM ) 

ENECCM-3) 

ENP(CM-3l 

XNUE(S-l) 

XNUI (S*l) 

0,00 

7.55E-06 

b.alE+03 

1 .53EM1 

6.00E+09 

5.00 

1 .83E-05 

7.66E*03 

8.22E  +  1  0 

3.39E+09 

10.00 

a.69£-05 

7.25E+03 

a,  11EM0 

l  .80E  +  09 

15,00 

5,62E-0a 

1 ,99£f0a 

1 .91E  +  1 0 

8.62E*08 

20.00 

1.57E-02 

8.24E+04 

8.72E+Q9 

3.93E*08 

25.00 

1 .57E-01 

1 .51E  +  05 

4.03E+09 

1 .80E  +  08 

30.00 

7 , 2 1 E*0 1 

1 ,  6  1 E  ♦  0  5 

1 .90E+09 

8 , 36E  +  07 

35.00 

3.13E+00 

1 .33E*05 

9,l9Et08 

3.95Et07 

ao.oo 

2.43e*01 

1  .00E  +  05 

4.62E+08 

1  .94E  +  07 

as. oo 

7.77E+01 

7 . 39E*0a 

2,aiEt08 

9.87E*06 

50.00 

1 .60E  +  02 

5.00E+04 

1 ,29Et08 

5.23ET06 

55.00 

8.00E+02 

2.23E+04 

6 , 96E+07 

3.26E+06 

60.00 

1 ,69E4’03 

7.93E+03 

3.65E+07 

2 , 0  1 E  +  06 

65.00 

2,a8E*03 

3.46£*03 

1 .86E  +  07 

1 ,2lEt06 

70.00 

2 . 5 1 E  +  0  3 

2.62E+03 

9.02E+06 

7.01E+05 

75.00 

2.«9e+03 

2,a9Ef 03 

4, 12E  +  06 

3.84E+05 

60,00 

3.07E+03 

3.07E+03 

1 ,75E*06 

1.96E+05 

85,00 

8.07e+03 

8.07Ef03 

6,99E*05 

9.21E+04 

90.00 

1 .36E*oa 

1 ,36E*0a 

2.73E+05 

4. 13E  +  04 

95.00 

3.77Etoa 

3.77etoa 

1 .10E  +  05 

1 ,87E*04 

100.00 

1 .00E  +  05 

1  ,OQE*05 

4.65E+04 

8.95E*03 

105.00 

1 . 95e ♦ 05 

1 .95E+05 

2, 14E+04 

3.95E+03 

110.00 

2.40E+05 

2.40ET05 

1 . 10E  +  04 

1.92E+03 

115,00 

2.20E+05 

2.20E405 

6.22E+03 

1 ,02E*03 

120.00 

2.a0E+05 

2.40E+05 

3.87E+03 

5.89E*02 

130.00 

3 . 06e+05 

3,06£*05 

1 .86E  +  03 

2,48Et02 

lao.oo 

2.a8E+05 

2,a«Ef 05 

1 .06E  +  03 

1 .33E  +  02 

160.00 

2 , 98E  t05 

2 , 98E ♦ 05 

4.32E+02 

5,27Et01 

180.00 

4.80E+05 

4.80E+05 

2.04E+02 

2.57E+0 1 

200.00 

1 ,09e*07 

1 .09E  +  07 

1 .04E  +  02 

1 ,4QE«-0  1 

Table  9.  Comparison  of  attenuation  rate  and  phase  velocity  from 
NOSC  and  Booker  models,  75  Hz. 


Booker  Model 

Profile  a  c/v 

NOSC  Model 
a  c/v 

WEDCOM 

Night 

2.1 

1.28 

1.5 

1.26 

NOSC 

Night 

1.8 

1.21 

0.9 

1.17 

WEDCOM 

Day 

1.6 

1.24 

1.7 

1.19 

NOSC 

Day 

1.3 

1.22 

1.3 

1.20 

Nucl ear 

4.3 

1.46 

3.8 

1.39 

Table  10.  Comparison  of  attenuation  rate  and  phase  velocity  from 
NOSC  and  Booker  models,  150  Hz. 


Booker  Model 

Profile  a  c/v 

NOSC  Model 
a  c/v 

WEDCOM 

Night 

3.1 

1.24 

2.5 

1.24 

NOSC 

Ni  ght 

2.9 

1.21 

1.3 

1.18 

WEDCOM 

Day 

2.7 

1.20 

2.2 

1.16 

NOSC 

Day 

2.5 

1.19 

2.2 

1.16 

Nuclear 

6.7 

1.34 

6.0 

1.29 
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where 


where  1^  is  an  effective  reflection  height.  Use  of  the  real  part  of  Zg  gives  good 
results  and  is  consistent  with  approximations  used  in  Reference  22. 

Tables  11  and  12  show  comparisons  of  height  gain  and  excitation  factor 
calculations  obtained  from  the  above  equations  (approximate  model)  and  the  full 
wave  model  described  in  Reference  14. 


Table  11.  Comparison  of  excitation  and  height  gain  factors  from  full 
wave  (first  entry)  and  approximate  (second  entry)  models. 


150  Hz. 


Profile 

ESI 

|ae| 

Jab  i 

leVl 

leEl 

m 

WEDCOM 

1.3 

1.6 

0.1 

2 

8 . 2  ( -  5 ) 

6. 2(-5) 

Ni  ght 

1.3 

1.7 

0 

2 

8 . 2  ( -  5 ) 

6 . 3  ( -  5 ) 

NO  SC 

1.2 

1.4 

0.1 

2 

8 . 2( - 5) 

5 . 3  ( -  5 ) 

Night 

1.3 

1.6 

0 

2 

8 . 2  ( -  5 ) 

5.8(-5) 

WEDCOM 

1.6 

1.8 

0.2 

2 

8.2(-5) 

5. 2( -5) 

Day 

1.6 

1.9 

0 

2 

8.2(-5) 

5 . 7  ( -  5 ) 

NOSC 

1.5 

1.7 

0.2 

2 

8.2(-5) 

5 . 1 ( -5 ) 

Day 

1.5 

1.8 

0 

2 

8 . 2  ( -  5 ) 

5 . 5  ( -  5 ) 

Nucl  ear 

2.2 

2.9 

0.1 

2 

8.2(-5) 

7 . 5  ( -  5 ) 

2.3 

3.2 

0 

2 

8 . 2( -5 ) 

8. 2(-5) 

Table  12.  Comparison  of  excitation  and  height-gain  factors  from  full 
wave  model  (first  entry)  and  approximate  (second  entry) 
models,  75  Hz. 


Profile 

n 

Hi 

IB 

H 

Ie£l 

|eBl 

WEDCOM 

2.7 

3.4 

0.3 

2 

5.8(-5) 

4.7(-5) 

Night 

2.6 

3.4 

0 

2 

5 . 8(  -  5 ) 

4 . 9  ( -  5 ) 

NOSC 

2.4 

2.8 

0.6 

2 

5.8(-5) 

3 . 6(  -  5 ) 

Nigh* 

2.5 

3.1 

0 

2 

5.8(-5) 

4 . 2(-5) 

WEDCOM 

3.7 

3.9 

0.4 

2 

5 . 8(  -  5 ) 

4 . 1  ( -  5 ) 

Day 

3.4 

4.2 

0 

2 

5.8(-5) 

4 . 5  ( -  5 ) 

NOSC 

3.3 

3.9 

0.5 

2 

5.8(-5) 

4 . 1  ( -  5 ) 

Day 

3.1 

3.8 

0 

2 

5.8(-5) 

4 . 2(  -  5 ) 

Nuclear 

5.1 

7.2 

0.4 

2 

5.8(-5) 

6 . 2  ( -  5 ) 

5.5 

8.1 

0 

2 

5.8(-5) 

6 . 8( - 5 ) 

SECTION  6 

CODE  STRUCTURE 


Major  computational  elements  used  in  link  vulnerability  and  simulation 
codes  include: 

•  Environment 

•  Propagation 

•  Signal  processing 

•  Network  status. 

The  network  status  calculation  can  be  included  in  simulation  codes  where  circuit 
performance  for  many  circuits  is  evaluated  or  performed  separately  using  results 
from  simulation  codes. 

A  desirable  goal  is  to  have  common  environment  calculations  for  communica¬ 
tion  codes.  Environment  includes  descriptions  of  the  natural  and  disturbed  atmos¬ 
pheres.  This  includes  electron  and  ion  densities  and  collision  frequencies  and  can 
include  dust,  water  vapor,  other  molecular  species,  and  possibly  terrain  modifica¬ 
tions  (cratering).  Link  vulnerability  codes  generally  use  engineering  models  for 
environment  calculations.  These  models  describe  specific  physical  processes  but 
use  engineering  approximations  for  part  or  all  of  the  phenomena.  The  models  require 
much  less  computational  time  than  first  principle  models  (eg,  detailed  MHD  calcula¬ 
tions)  but  still  provide  sufficient  detail  and  accuracy  for  systems  applications. 
Multiburst,  multilink  simulation  codes  require  even  faster  running  models  in  order 
to  handle  the  large  spatial  and  temporal  regions  required.  This  is  accomplished 
through  the  use  of  simplified  algorithms  or  a  precomputed  data  base  that  can  be 
interpolated.  By  developing  the  simulation  environment  models  from  those  developed 
for  link  vulnerability  codes  the  desired  consistency  between  the  communication  code 
environment  models  can  be  maintained. 

The  engineering  ionization  calculations  for  link  vulnerability  codes  can 
be  conveniently  divided  into  the  following  models: 

•  Natural  (undisturbed)  atmosphere 

•  Natural  ionosphere 

•  Fireball/debris  geometry 

•  D-region  ionization 


•  E-  and  F-region  ionization 

•  Fireball  ionization. 

In  performing  these  calculations  the  spatial  region  requirements  and  the 
ionization  modeling  are  significantly  different  for  ionospheric  dependent  propaga¬ 
tion  (ELF  through  HF)  and  1 ine-of-sight  (LOS)  propagation  in  and  above  the  VHF  band. 
In  the  latter  case  ionization  calculations  are  needed  along  LOS  paths.  In  general 
the  LOS  paths  can  move  with  time  and  calculations  must  either  be  made  along  each 
path  of  interest  or  a  three-dimensional  grid  and  interpolation  procedure  must  be 
used.  If  D-region  ionization  is  to  be  included,  the  grid  size  must  be  relatively 
small  or  calculations  are  required  along  the  LOS  path  in  the  D  region  and  the  grid 
only  used  in  the  E  and  F  regions  (method  used  in  ROSCOE) .  Generally,  equilibrium 
ionization  calculations  can  be  used  for  ionization  levels  affecting  propagation  in 
and  above  the  VHF  band.  Finally,  it  is  generally  efficient  to  combine  ionization 
and  propagation  calculations  for  LOS  paths. 

For  ionospheric  dependent  propagation,  calculations  are  required  in  a 
vertical  plane  defined  by  the  great  circle  path  between  transmitter  and  receiver 
(this  neglects  off-great  circle  paths).  For  ELF  and  VLF  propagation,  ionization 
along  vertical  profiles  is  required  for  mode  solutions.  For  LF  and  HF  propagation 
ionization  calculations  are  required  along  skywave  paths.  Since  a  relatively  large 
number  of  skywaves  can  propagate,  it  is  generally  efficient  to  determine  the  ioniza¬ 
tion  along  vertical  profiles  and  then  use  these  to  determine  ionization  along  the 
skywaves.  The  ionization  levels  that  can  affect  ionospheric  dependent  propagation 
are  small  enough  that  nonequilibrium  conditions  can  be  important.  These  include 
the  buildup  of  ionization  after  a  debris  region  has  moved  into  an  area  and  the  decay 
of  ionization  after  a  debris  region  has  moved  away.  For  ionospheric  dependent  prop¬ 
agation  it  is  generally  not  efficient  to  combine  ionization  and  propagation 
calculations . 

For  both  LOS  propagation  and  ionospheric  dependent  propagation  the  fireball 
and  debris  geometry  can  be  calculated  from  the  same  model.  For  LOS  propagation  fire¬ 
ball  ionization  along  those  paths  that  intersect  the  fireball  can  be  computed  from 
fireball  geometry  and  ionization  models.  For  ionospheric  dependent  propagation 
fireball  ionization  along  the  vertical  paths  used  to  define  the  plane  between  ter¬ 
minals  can  be  determined  but  ionization  affecting  off-great  circle  paths  is  difficult 
to  determine  in  an  efficient  manner.  Use  of  simplified  spatial  distributions  may  be 
adequate  for  estimates  or  propagation  flags. 
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Figure  34  shows  simplified  flow  diagrams  for  the  ionization  models  re¬ 
quired  for  LOS  and  ionospheric  dependent  propagation  calculations.  The  time  loop 
could  be  placed  outside  the  path  loop  for  ionospheric  dependent  propagation  but 
this  requires  considerable  storage  in  order  to  perform  nonequilibrium  ionization 
calculations . 

The  selection  of  the  vertical  path  locations  for  ionospheric  dependent 
ionization  calculations  will  affect  both  accuracy  and  computation  time.  For  most 
natural  ionospheric  conditions  vertical  path  spacings  of  several  thousands  of 
kilometers  are  adequate.  Smaller  spacings  are  generally  required  for  disturbed 
conditions.  A  reasonably  efficient  method  for  determining  the  vertical  path  loca¬ 
tions  appears  to  be  to  first  select  locations  based  on  the  natural  ionosphere. 

Then  for  each  interval  between  vertical  paths  tests  could  be  made  to  see  if  nuclear 
ionization  is  important.  If  it  is,  the  interval  could  be  divided  into  smaller 
intervals.  In  order  to  determine  whether  nuclear  ionization  is  important  the  great 
circle  distance  from  each  terminal  and  from  the  center  of  the  interval  between 
terminals  to  the  burst  point  (prompt  ionization)  or  to  the  debris  center  (delayed 
ionization)  would  be  determined  and  compared  to  an  effects  radius.  The  effects 
radius  would  be  a  function  of  burst  altitude  (prompt  ionization)  or  debris  altitude 
(delayed  ionization)  and  frequency.  Since  the  nuclear  ionization  can  move  with 
time,  the  tests  would  be  made  at  all  calculation  times  and  an  interval  spacing 
between  vertical  paths  would  only  remain  at  the  value  chosen  for  the  natural  ionos¬ 
phere  if  nuclear  ionization  was  not  important  at  any  of  the  calculation  times.  An 
interval  spacing  of  2000  km  for  natural  conditions  and  500  km  for  disturbed  condi¬ 
tions  appears  consistent  with  accuracy  and  computational  time  requirements.  A 
smaller  spacing  for  regions  affected  by  surface  or  near-surface  bursts  may  be  re¬ 
quired  for  some  systems  applications  but  generally  ionospheric  dependent  propagation 
is  not  significantly  affected  by  small  ionization  regions. 

The  altitude  interval  used  to  compute  ionization  for  the  vertical  paths 
depends  on  propagation  frequency.  For  ELF  propagation,  ionization  calculations  are 
required  below  about  200  km.  For  VLF  and  LF  propagation  ionization  calculations  are 
required  below  about  100  km.  For  HF  propagation  the  altitude  interval  required 
depends  on  whether  skywave  geometries  are  determined  from  precomputed  results  for 
natural  conditions  or  are  found  by  ray  tracing.  If  precomputed  results  are  used, 
ionization  calculations  are  only  needed  below  about  120  km.  Ray  tracing  would  re¬ 
quire  ionization  calculations  below  about  400  km. 


Several  methods  are  used  to  select  calculation  times  in  communication 
codes.  In  the  current  WESCOM  and  WEDCOM  codes  the  calculation  times  are  selected 
from  user  input  data.  In  HFNET  calculations  are  made  at  update  times  determined  by 
the  code  and  then  calculations  at  event  times  (user  requested  times  or  times  when 
the  codes  determines  that  calculations  are  required)  are  found  by  interpolation. 

The  use  of  user  selected  calculation  times  is  the  simplest  method  of  time  selection 
since  it  does  not  require  update  times  or  interpolation  procedures.  However,  it 
requires  that  the  user  have  some  knowledge  of  significant  times  and  reasonable  time 
intervals  between  calculations.  The  use  of  code  selected  times  allows  determining 
times  based  on  environment  calculations.  If  an  ionization  file  is  prepared  for  code 
selected  update  times,  it  can  be  subsequently  used  to  determine  ionization  at  user 
selected  calculation  times.  By  saving  the  ionization  file  the  user  selected  times 
could  be  changed  in  subsequent  calculations. 

A  three-level  option  for  selection  of  calculation  times  appears  useful. 

The  first  option  would  be  the  use  of  user  selected  times  as  in  the  current  WEDCOM 
code.  The  second  option  would  be  the  use  of  code  determined  times  between  input 
start  and  stop  times.  The  third  option  would  be  the  use  of  the  ionization  file 
prepared  for  option  2  to  determine  ionization  at  user  selected  times. 
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S.  Gutsche 
R.  Hendrick 

G.  McCartor 
Tech  Library 
F.  Guigliano 


Mitre  Corp 

ATTN:  A.  Kymmel 
ATTN:  C.  Callahan 
ATTN  MS  J104 ,  M.  Dresp 
ATTN:  G.  Harding 
ATTN:  B.  Adams 
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DEPARTMENT  OF  DEFENSE  CONTRACTORS  (Continued) 


DEPARTMENT  OF  DEFENSE  CONTRACTORS  (Continued) 


Mitre  Corp 

Science  Applications,  Inc 

ATTN: 

M.  Horrocks 

ATTN 

SZ 

ATTN: 

W.  Hall 

ATTN: 

J.  Wheeler 

Science  Appl  ications,  Inc 

ATTN: 

W.  Foster 

ATTN 

J.  Cockayne 

Pacific-Sierra  Research  Corp 

SRI  International 

ATTN: 

F.  Thomas 

ATTN 

G.  Price 

ATTN: 

E.  Field,  Jr 

ATTN 

R.  Tsunoda 

ATTN: 

H.  Brode,  Chairman  SAGE 

ATTN 

J.  Vickrey 

ATTN 

W.  Chesnut 

Pennsyl vania 

State  University 

ATTN 

D.  Neil  son 

ATTN: 

Ionospheric  Rsch  Lab 

ATTN 

J.  Petrickes 

ATTN 

R.  Leada brand 

Photometries 

Inc 

ATTN 

R.  Livingston 

ATTN: 

I.  Kofsky 

ATTN 

M.  Baron 

ATTN 

A.  Burns 

Physical  Dynamics,  Inc 

ATTN 

C.  Rino 

ATTN: 

E.  Fremouw 

ATTN 

W.  Jaye 

Physical  Rsch,  Inc 


ATTN: 

R. 

Del iberis 

Associates 

ATTN: 

W. 

Karzas 

ATTN: 

W. 

Wright 

ATTN: 

C. 

Greif inger 

ATTN: 

R. 

Lelevier 

ATTN: 

F. 

Gilmore 

ATTN: 

B. 

Gabbard 

ATTN: 

H. 

Ory 

ATTN: 

R. 

Turco 

ATTN: 

M. 

Gantsweg 

Assoc 

iates 

ATTN: 

B. 

Yoon 

Rand  Corp 

ATTN:  E.  Bedrozian 
ATTN:  C.  Crain 

Riverside  Rsch  Institute 
ATTN:  V.  Trapani 

Rockwell  International  Corp 
ATTN:  S.  Qullicl 

Santa  Fe  Corp 

ATTN:  D.  Paoluccl 

Science  Applications,  Inc 


Rockwell  International  Corp 
ATTN:  R.  Buckner 


ATTN 

ATTN 

ATTN 


G.  Smith 
V.  Gonzales 
D.  McDaniels 


ATTN: 

L.  Linson 

ATTN: 

J.  Carpenter 

ATTN: 

E.  Straker 

ATTN: 

C.  Humphrey 

ATTN: 

D.  Hamlin 

ATTN: 

W,  Reldy 

ATTN: 

C.  Smith 

ATTN: 

0.  Shepard 

Stewart  Radiance  Lab 
ATTN:  J.  Ulwich 

Syl vania  Systems  Group 

ATTN:  R.  Steinhoff 

Strategic  Systems  Dlv 

ATTN:  J.  Concordia 
ATTN:  I.  Kohl  berg 

Technology  International  Corp 
ATTN:  H.  Boquist 

Tri-Corn,  Inc 

ATTN:  D.  Murray 

TRW  Electronics  S  Defense  Sector 
ATTN:  D.  Dee 
ATTN:  R.  Plebuch 

Utah  State  University 

ATTN:  Sec  Con  Ofc  for  D.  Burt 

ATTN:  Sec  Con  Ofc  for  K.  Baker, 

Dir  Atmos  A  Space  Scl 

ATTN:  Sec  Con  Ofc  for  L.  Jensen,  Elec  Eng  Dept 

ATTN:  Sec  Con  Ofc  for  a.  Steed 

Vlsldyne,  Inc 
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